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Abstract 

The  reliability  of  a  single-unit  system  experiencing  degradation  (wear)  due  to 
the  influence  of  a  general,  observable  environment  process  is  considered.  In  partic¬ 
ular,  the  failure  time  distribution  is  evaluated  using  only  observations  of  the  unit’s 
current  operating  environment  which  is  characterized  as  a  finite  semi-Markov  process 
(SMP).  In  order  to  impose  the  Markov  property,  generally  distributed  environment 
state  sojourn  times  are  approximated  as  phase-type  (PH)  random  variables  using 
observations  of  state  holding  times  and  transition  rates.  The  use  of  PH  distributions 
facilitates  the  use  of  existing  analytical  results  for  reliability  evaluation  of  units  sub¬ 
ject  to  an  environment  process  that  evolves  as  a  continuous-time  Markov  chain.  The 
procedure  is  illustrated  through  three  numerical  examples,  and  results  are  compared 
with  those  obtained  via  Monte  Carlo  simulation.  The  maximum  absolute  deviation 
in  probability  for  failure  time  distributions  was  on  the  order  of  0.004.  The  results  of 
this  thesis  provide  a  novel  approach  to  the  reliability  analysis  of  units  operating  in 
randomly  evolving  environments  for  which  degradation  or  failure  time  observations 
are  difficult  or  impossible  to  obtain. 
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PHASE-TYPE  APPROXIMATIONS  FOR  WEAR  PROCESSES  IN 
A  SEMI-MARKOV  ENVIRONMENT 


1.  Introduction 

1 . 1  Background, 

Classical  reliability  analysis  involves  the  estimation  of  the  probability  distribu¬ 
tion  of  component  lifetimes  based  on  historical  observations  of  failure  time.  This  has 
been  the  most  thoroughly  studied  form  of  reliability  analysis,  and  several  standard 
models  have  been  used  for  decades  in  the  design,  manufacture,  and  maintenance  of 
components  and  systems.  As  an  example,  the  exponential  probability  distribution 
has  been  widely  accepted  as  an  estimate  for  the  lifetime  distribution  of  electronic 
components.  However,  advances  in  technology  and  production  methods  have  made 
it  impractical  to  observe  failures  of  systems  with  very  long  lifetimes.  Additionally, 
the  observation  of  failure  times  may  be  a  prohibitively  expensive  approach  given 
the  complexity  and  high  costs  of  modern  systems.  For  instance,  consider  a  satellite 
operating  in  orbit.  Evaluating  the  reliability  of  a  single  component  or  the  entire 
system  by  running  until  failure  is  clearly  an  infeasible  approach.  An  alternative 
technique  commonly  applied  to  alleviate  these  shortfalls  in  failure  time  estimation  is 
accelerated  lifetime  testing.  This  technique  involves  lifetime  testing  of  components 
under  accelerated  laboratory  conditions.  However,  the  value  of  this  approach  may 
be  countered  by  the  fact  that  such  artificial  laboratory  tests  often  do  not  accurately 
represent  the  realistic  conditions  of  a  unit’s  operating  environment. 

Extremely  long  lifetimes,  high  costs,  and  unrealistic  testing  environments  are 
not  the  only  drawbacks  of  failure-based  reliability  analysis.  A  system’s  reliability 
may  not  necessarily  be  defined  as  the  time  of  a  catastrophic  failure-achievement  of  a 
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threshold  level  of  wear  or  degradation  may  signal  the  need  for  maintenance  or  the  end 
of  a  system’s  useful  operating  life.  With  failure-based  reliability  analysis,  these  mea¬ 
sures  are  not  necessarily  observed.  In  an  effort  to  overcome  these  shortfalls,  a  more 
recent  trend  in  system  lifetime  estimation  has  been  the  application  of  degradation- 
based  reliability  techniques.  In  contrast  to  failure-based  approaches,  these  techniques 
involve  lifetime  estimation  of  a  component  or  system  based  on  observations  of  cumu¬ 
lative  degradation  over  time.  However,  degradation  measurements  may  be  difficult 
or  even  impossible  to  obtain.  For  instance,  it  may  be  infeasible  or  prohibitively 
expensive  to  record  degradation  measurements  for  a  solar  panel  attached  to  an  on- 
orbit  satellite.  Although  this  technique  may  more  accurately  predict  the  degradation 
path  leading  to  system  failure,  it  does  not  necessarily  take  into  account  the  ambient 
environment  that  is  causing  the  unit  to  degrade. 

In  such  cases  as  the  on-orbit  satellite,  the  environment  to  which  the  system  is 
exposed  is  of  great  importance.  For  instance,  differing  levels  of  ultraviolet  radiation 
may  have  varied  effects  on  the  overall  wear  of  the  solar  panel.  An  environment- 
based  approach  to  reliability  analysis  seeks  to  account  for  the  impact  the  surrounding 
environment  has  on  the  operation  and  degradation  of  the  system  under  study.  Since 
very  few  systems  operate  in  complete  isolation,  environment-based  approaches  are 
valuable  in  evaluating  reliability  under  realistic  operating  conditions.  In  lieu  of 
imposing  an  artificial  setting,  environment-based  reliability  models  attempt  to  study 
both  the  operating  environment  and  its  effects  on  the  system.  As  with  accelerated 
testing  techniques,  this  approach’s  success  in  estimating  system  reliability  depends  on 
the  level  of  accuracy  of  the  environment  model  under  which  the  system  is  studied  and 
the  effect  of  various  environment  conditions  on  the  degradation  of  the  system.  If  the 
environment  consists  of  easily  recognizable  and  definable  discrete  states,  modelling 
of  the  environment  may  be  a  fairly  straightforward  process.  However,  characterizing 
and  defining  the  various  states  of  a  continuous  environment  process  may  not  be  an 
easy  task.  Furthermore,  in  some  cases,  the  transitions  from  one  environment  state 
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to  another  are  not  easily  observable.  Therefore,  another  potential  drawback  of  an 
environment-based  reliability  approach  is  in  the  discretization  of  its  feasible  state 
space.  Finally,  it  may  be  difficult  or  even  impossible  to  observe  the  environment. 
Systems  operating  under  deep  water  or  in  outer  space  are  examples  in  which  the 
environment  may  be  difficult  to  characterize. 

Despite  the  drawbacks  noted  above,  environment-based  reliability  models  have 
potential  for  accurate  lifetime  estimation  of  systems  operating  in  known  or  observ¬ 
able  environments.  Several  techniques,  both  probabilistic  and  statistical  in  nature, 
have  been  developed  to  evaluate  the  reliability,  and  particularly  the  failure  time 
distribution,  of  a  system  experiencing  wear  due  to  a  dynamic  environment  process. 
According  to  Singpurwalla  [25],  an  advantage  of  probabilistic  approaches,  via  the 
use  of  stochastic  models,  is  that  both  the  physics  of  system  or  component  failure 
and  the  dynamic  nature  of  the  environment  are  taken  into  account.  Following  this 
approach,  Kharoufeh  [16]  studied  the  reliability  (via  derivation  of  the  failure  time 
distribution)  of  a  single-unit  system  experiencing  wear  due  to  a  multi-state  stochas¬ 
tic  environment  process.  This  technique  uses  a  Markov  additive  process  (MAP),  as 
reviewed  by  Singpuwalla  [25].  Using  the  MAP  as  the  basis  for  the  overall  model, 
the  cumulative  damage  of  the  system  is  modelled  as  a  continuous,  nondecreasing 
stochastic  process,  and  the  random  environment  process  is  modelled  as  a  finite  state 
space  stochastic  process.  In  [16],  Kharoufeh  assumed  the  environment  process  to  be 
a  continuous-time  Markov  chain  (CTMC)  in  which  the  evolution  of  the  environment 
depends  only  upon  the  current  state,  and  the  time  spent  in  each  state  is  assumed  to 
be  an  exponentially  distributed  random  variable,  ffowever,  this  model’s  assumption 
of  exponentially  distributed  state  holding  times  limits  its  real-world  applicability. 
In  many  situations,  the  distribution  of  the  random  time  spent  in  each  state  of  the 
environment  may  be  either  nonexponential  or  simply  unknown.  A  probabilistic  tech¬ 
nique  based  solely  on  measurements  of  the  environment  state  holding  times  would 
alleviate  the  restrictive  need  for  exponential  (or  any  other)  state  holding  time  dis- 
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tributions.  This  thesis  effort,  based  on  the  environment-based  approach,  attempts 
to  address  this  shortcoming  by  presenting  numerical  approximation  procedures  for 
general  environment  models. 

1.2  Problem  Definition  and  Methodology 

This  thesis  proposes  a  method  for  estimating  the  lifetime  of  a  system  based 
solely  on  observations  of  the  environment  to  which  the  system  is  exposed  and  knowl¬ 
edge  of  the  environment’s  impact  on  the  degradation  of  the  system.  In  particular, 
the  reliability  of  a  single-unit  system  experiencing  nondecreasing  wear  due  to  a  multi¬ 
state  environment  process  in  which  the  evolution  of  the  environment  process  depends 
only  upon  its  current  state  is  investigated.  Furthermore,  it  is  assumed  that  the  time 
spent  in  each  state  is  known  only  to  be  a  positive-valued,  nondefective  random  vari¬ 
able.  That  is,  there  is  zero  probability  that  the  environment  process  will  remain  in 
any  one  state  indefinitely. 

This  work  will  build  upon  and  generalize  the  results  of  [16]  and  [17].  In  those 
works,  the  authors  assumed  the  wear-inducing  environment  to  be  a  continuous-time 
Markov  chain  (CTMC).  As  such,  the  next  state  of  the  environment  depends  ex¬ 
clusively  on  the  current  state,  forming  an  embedded  discrete-time  Markov  chain 
(DTMC).  Additionally,  the  time  spent  in  each  of  the  various  environment  states  is 
assumed  to  be  an  exponentially  distributed  random  variable.  In  reality,  state  hold¬ 
ing  time  distributions  may  not  follow  exponential  distributions  at  all.  In  fact,  the 
particular  state  holding  time  distributions  may  be  completely  unknown.  However, 
it  may  be  assumed  that  transitions  to  future  states  of  the  environment  are  still  de¬ 
pendent  only  upon  the  current  state,  thus  retaining  an  embedded  DTMC.  In  such 
cases,  the  environment  process  may  be  characterized  as  a  semi-Markov  process.  In 
this  work,  the  primary  focus  is  the  reliability  analysis  of  systems  experiencing  cumu¬ 
lative  wear  damage  due  to  the  influence  of  a  dynamic  external  environment  modelled 
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as  a  finite  state  space  semi-Markov  process  which  allows  for  general  state  holding 
time  distributions. 

The  case  of  nonexponential  environment  state  holding  times  is  studied  by  ap¬ 
proximating  arbitrary  distributions  with  phase-type  (PH)  distributions.  That  is,  the 
arbitrary  holding  time  distributions  are  approximated  as  the  time  to  absorption  of 
an  underlying  Markov  process  that  must  be  estimated.  In  order  to  construct  such 
approximations,  the  environment  process  must  be  observed  for  some  period  of  time, 
and  the  time  spent  in  each  environment  state  must  be  recorded.  Furthermore,  the 
number  of  transitions  among  the  various  states  must  be  observed  and  recorded  in 
order  to  statistically  estimate  the  inter-state  transition  rates.  The  advantage  of  PH 
distributions  is  that  they  retain  the  Markov  property,  supplanting  the  need  for  ex¬ 
ponential  state  holding  time  distributions.  This  thesis  will  demonstrate  that,  by 
imposing  the  Markovian  property  by  means  of  phase-type  approximations,  the  semi- 
Markov  environment  process  can  be  converted  to  a  CTMC,  and  the  results  developed 
by  Kharoufeh  [16]  and  Kharoufeh  and  Sipe  [17]  can  be  applied. 

The  results  of  this  research  can  be  applied  to  the  reliability  analysis  of  com¬ 
ponents  or  systems  that  degrade  over  time  due  to  exposure  to  varying  levels  of 
temperature,  pressure,  or  some  other  environmental  measure.  This  thesis  effort  will 
be  valuable  to  those  who  design,  evaluate,  repair,  or  replace  systems  exposed  to  ran¬ 
dom  environments.  In  particular,  designers  of  components  used  aboard  satellites  or 
underwater  vessels,  systems  known  to  operate  in  extreme  environments,  may  find  the 
technique  developed  in  this  thesis  to  be  particularly  useful.  This  work  will  also  ad¬ 
vance  the  current  state-of-the-art  in  estimation  techniques  via  stochastic  modelling 
of  degradation  and  environment  processes. 
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1.3  Outline  of  the  Thesis 


The  remainder  of  this  thesis  is  organized  as  follows.  Chapter  2  reviews  the  cur¬ 
rent  literature  on  several  probabilistic  and  statistical  degradation-based  reliability 
analysis  techniques.  Chapter  3  describes  the  formal  mathematical  model,  presents 
the  techniques  used  to  construct  phase-type  approximations  of  arbitrary  distribu¬ 
tions,  and  proposes  a  method  by  which  single-unit  lifetime  distributions  may  be  ap¬ 
proximated  in  semi- Markovian  environments.  Chapter  4  compares  the  failure  time 
distribution  results  of  numerical  examples  of  the  proposed  technique  with  those  ob¬ 
tained  via  simulation.  Finally,  chapter  5  presents  conclusions  and  some  suggestions 
for  future  research. 
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2.  Review  of  the  Literature 

This  chapter  provides  an  overview  of  the  literature  concerning  techniques  for 
the  analysis  of  degradation-  or  environment-based  reliability.  Although  certainly  not 
an  exhaustive  survey,  it  attempts  to  provide  a  sense  of  the  two  prominent  methods 
for  evaluating  the  lifetime  of  systems  experiencing  wear  over  time:  probabilistic  and 
statistical  approaches. 

2.1  Probabilistic  Approaches  to  Degradation- Based  Reliability 

Several  probabilistic  approaches  to  degradation-based  reliability  are  found  in 
the  literature.  In  his  extensive  survey  of  probabilistic  approaches  to  failure  time 
analysis,  Singpurwalla  [25]  considers  both  the  physics  of  failure  and  the  characteris¬ 
tics  of  the  operating  environment.  In  general,  the  variation  of  the  degradation  rate  is 
caused  by  a  dynamic  environment;  this  environment  can  be  described  by  a  stochastic 
process.  However,  failure  time  distributions  derived  from  stochastic  processes  often 
do  not  have  closed-form  expressions  and  can  only  be  expressed  via  Laplace  trans¬ 
forms,  as  seen  in  [16].  Singpurwalla  [25]  notes  four  broad  strategies  that  have  been 
developed  for  modelling  failure  distributions  based  on  stochastic  processes: 

1.  Wear  is  described  by  a  stochastic  process  with  continuous  sample  path  (i.e.  a 
diffusion  process),  such  as  Wiener,  gamma,  or  deterministic; 

2.  failure  rate  is  described  by  a  stochastic  process,  such  as  gamma,  shot-noise, 
functionals  of  a  Wiener  process,  or  a  Levy  process; 

3.  the  environment  is  described  by  a  stochastic  process,  typically  a  shock-inflicting 
Poisson  process; 

4.  a  response  variable  correlated  with  the  lifclcngth  is  described  by  a  stochastic 
process,  such  as  a  stationary,  continuous-time  Gaussian  process. 
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When  wear  is  modelled  as  a  stochastic  process  with  continuous  sample  path,  Singpur- 
walla  [25]  notes  that  use  of  a  Wiener  process  to  describe  diffusion  of  the  item  state 
may  provide  better  goodness-of-fit  than  use  of  the  exponential  and  Weibull  distribu¬ 
tions.  In  [25],  the  author  describes  such  a  process,  also  known  as  Brownian  motion, 
as  follows.  The  process  (7(5),  s  >  0}  is  such  that  7(0)  =  0,7(5)  has  stationary 
independent  incrememcnts,  and,  for  every  s  >  0,  7(5)  ~  A/"(0,c2s)  for  some  c  >  0. 
Singpurwalla  [25]  adds  that  the  Wiener  process  is  a  Markov  process  since  the  inde¬ 
pendent  increments  property  of  {7(5),  s  >  0}  implies  that,  for  any  t  >  0, 

P{l(t  +  s)  <  a  j  7 (s)  =  x ,  7 (u),  0  <  u  <  s}  =  P{7(t  +  s)  <  a  |  7(5)  =  x}.  (2.1) 

According  to  Singpurwalla  [25],  a  disadvantage  of  the  use  of  a  Wiener  process  to 
model  the  diffusion  of  item  state  is  that  the  sample  path  of  wear  is  not  monotonically 
increasing. 

As  an  example  of  the  failure  rate  (i.e.  hazard  rate  function)  being  described 
by  a  stochastic  process,  Singpurwalla  [25]  considers  a  shot  noise  hazard  rate  process 
as  follows.  Suppose  an  item  operates  in  a  dynamic  environment  that  inflicts  shocks 
over  time  u  according  to  a  Poisson  process  with  rate  m(u),  u  >  0.  Each  shock  results 
in  a  stress  on  the  item  that  contributes  to  its  failure  rate.  If  a  shock  of  magnitude 
D  occurs  at  an  epoch  S,  and  h(t),  a  nonincreasing  function  of  t,  is  the  attenuation 
function,  then  the  contribution  of  the  shock  to  the  item’s  failure  rate  is  Dh(t )  at 
time  S  + 1.  The  author  presents  the  item’s  failure  rate  at  time  t  as 

OO 

A  (t)  =  J2Dnh(t-Sn),  (2.2) 

n=l 

where  {Sn,  n  >  1}  are  the  epochs  at  which  the  shocks  occur,  {Dn,n  >  1}  are  the 
respective  shock  magnitudes,  and  h{t)  =  0  when  t  <  0. 

In  the  case  where  the  wear-inducing  environment  is  modelled  as  a  stochastic 
process,  Singpurwalla  [25]  summarizes  the  use  of  Markov  additive  processes  (MAP). 
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As  an  example  of  a  MAP,  one  studies  two  stochastic  processes:  X(t),  which  repre¬ 
sents  the  cumulative  damage  of  the  item  at  time  t,  and  Z(t ),  which  represents  the 
state  of  the  environment  at  time  t.  This  MAP  approach  was  used  by  Kharoufeh  [16] 
and  Kharoufeh  and  Sipe  [17]  and  is  reviewed  in  detail  in  the  next  chapter.  Addition¬ 
ally,  a  wear-inducing  shock  process  may  be  incorporated  into  this  model.  Singpur- 
walla  [25]  suggests  that  the  MAP  modelling  approach,  along  with  other  stochastic 
modelling  techniques,  will  become  most  useful  only  upon  the  further  development  of 
computational  methodologies  for  reliability  and  survival  analysis. 

In  the  fourth  and  final  strategy  for  the  development  of  failure  models  studied 
in  [25],  the  author  briefly  discusses  the  approach  in  which  a  response  variable  cor¬ 
related  with  the  lifelength  is  described  by  a  stochastic  process,  such  as  a  stationary, 
continuous-time  Gaussian  process.  According  to  the  author,  there  has  been  little  de¬ 
velopment  in  this  area  of  research.  Overall,  Singpurwalla  [25]  emphasizes  the  value 
of  stochastic-based  approaches  for  developing  failure  models,  concluding  that  these 
models  better  exploit  the  physics  of  failure  process.  One  drawback,  as  noted  in  [25], 
is  the  complexity  of  the  resulting  distributional  forms,  which  must  often  be  presented 
via  Laplace  transforms. 

Abdcl-Hameed  [2]  studies  properties  of  the  failure  time  distribution  when  the 
wear  process  is  assimied  to  be  an  increasing  Levy  process.  Qinlar  [10]  introduced 
variations  of  a  general  shock  and  wear  model  where  shocks  of  random  magnitude 
occur  at  environment  state  transitions  and  cause  degradation  to  the  system.  In 
[10],  the  author  incorporates  the  wear  and  environment  processes  into  a  MAP  and 
includes  examples  in  which  the  wear  process  is  both  a  Gamma  and  an  increasing 
Levy  process.  Esary  et  al.  [12]  consider  the  failure  time  distribution  of  a  system 
subject  to  shocks  governed  by  a  Poisson  process  and  then  generalize  the  problem  for 
cases  of  continuous  wear.  Linmios  and  Oprisan  [19]  apply  Markov  renewal  theory  to 
semi-Markov  systems  to  obtain  reliability  and  availability  measures. 
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Kharoufeh  [16]  provides  a  brief  overview  of  probabilistic  techniques  and  notes 
that  stochastic  shock  and  wear  models  are  the  most  commonly  studied  approaches. 
In  the  reliability  analysis  of  a  system  that  degrades  over  time,  the  author  character¬ 
izes  the  continuous  wear  process  as  a  nondecreasing  stochastic  process.  In  that  work, 
the  author  considered  a  single-unit  system  whose  cumulative  damage  over  time  is  a 
continuous  wear  process  dependent  on  an  external  environment  process  assumed  to 
be  a  temporally  homogeneous,  Markov  process  with  a  finite  state  space.  Kharoufeh 
[16]  derived  the  failure  time  distribution  and  moments  in  terms  of  Laplace-Stieltjes 
transforms.  Further  discussion  of  the  results  in  [16]  is  presented  in  the  next  chapter. 

2.2  Statistical  Approaches  to  Degradation- Based  Reliability 

In  addition  to  the  probabilistic  models  described  above,  the  literature  contains 
models  for  degradation-based  failure  distributions  derived  via  nonlinear  regression 
and  other  statistical  techniques.  Ahmad  and  Sheikh  [3]  present  characteristics  and 
applications  of  the  Bernstein  probability  density  function  (pdf).  Also  known  as  the 
tt-distribution,  or  inverted  normal  distribution,  the  Bernstein  pdf  was  first  developed 
by  Gertsbakh  and  Kordonsky  [14]  and  Vysokovskii  [27]  to  model  the  life  character¬ 
istics  of  machine  components  experiencing  non-stationary  linear  wear.  According  to 
Ahmad  and  Sheikh  [3],  the  Bernstein  probability  distribution  has  been  successfully 
used  to  model  cutting  tool  life,  monitoring  the  dimensions  of  machine  parts,  testing 
the  slideways  and  rotating  parts  of  machine  tools,  and  determining  tool  replacement 
intervals  in  precision  machining.  A  potential  drawback  of  the  Bernstein  probability 
distribution,  with  respect  to  its  use  in  reliability  analysis,  is  that  it  belongs  to  a 
bimodal  family  of  probability  distributions  for  which  no  moments  exist  [3].  This 
reliability  model  is  based  on  a  random  wear  process  of  the  form 

W(t)  =  at  +  b,  (2.3) 
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where  a  is  a  random  variable  that  describes  the  rate  of  wear,  W (t),  and  b  =  W(0)  is 
a  random  variable  denoting  the  initial  value  of  wear.  When  T  is  the  lifetime  of  the 
component,  and  W\  is  the  permissible  wear  limit,  the  three-parameter  Bernstein’s 
distribution  is  given  by 


F{t)  =  P(T  <  t)  =  $ 


t  —  c 

\J  oit2  +  (3 


(2.4) 


where 


<3>(1/a /a)  ~  1, 


-  E(b) 
C~  E(a)  ’ 

a  =  V(a)/E2(a ), 


(3  =  V(b)/E2(a), 


and 


c,  a,  (5,  E[(a)}  >  0. 

Ahmad  and  Sheik  [3]  further  propose  a  two-parameter  inverted  family  of  dis¬ 
tributions  with  (3  —  0  to  model  wear-related  life.  Their  work  includes  a  thorough 
study  of  the  properties  of  the  Bernstein  distribution,  including  maximum  likelihood 
estimates  of  its  parameters.  Lu  and  Meeker  [20]  expand  the  statistical  study  of 
degradation-based  time-to-failure  distributions  by  reviewing  methods  for  fitting  non¬ 
linear  regression  models  to  observed  fatigue-crack-growth  measurements.  The  pa¬ 
rameters  of  a  mixed-effect  path  model  are  estimated  using  a  two-stage  method.  The 
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authors  present  a  parametric  degradation  model  with  hxed-effect  parameters  com¬ 
mon  for  all  units  inspected  and  random-effect  parameters,  representing  individual 
unit  characteristics.  Since  expressing  the  time-to-failure  distribution  can  be  compli¬ 
cated  by  the  inclusion  of  more  than  one  random  parameter,  the  authors  proposed 
the  use  of  Monte  Carlo  simulation  to  obtain  the  distribution  numerically.  For  a 
degradation  path  modeled  by 


n(t)  =  0  +  0(j), 


(2.5) 


where  0  is  fixed  and  represents  the  common  initial  amount  of  degradation  of  all  test 
units,  0,  which  represents  the  degradation  rate,  varies  from  unit  to  unit  according 
to  a  Weibull(a,  (3)  distribution,  and  D  is  the  critical  degradation  level,  the  time-to- 
failure  distribution  is  given  by 


FT{t)  =  exp 


D-d 

at 


0 


,  t  >  0. 


(2.6) 


Since  the  random  variable  1/T  follows  a  Weibull  distribution,  this  time-to-failure 
distribution  is  known  as  the  reciprocal  Weibull.  When  0  follows  a  lognormal  distri¬ 
bution  with  parameters  (/i,<r2),  T  follows  a  lognormal  distribution,  and 


Frit) 


$ 


log(t)  -  \log(D  -  0)  ~  fA 
cr 


(2.7) 


When  0  ~  a2)  and  cr  -C  /x  (to  make  P{@  <  0}  negligible),  the  time-to-failure 

approximates  a  special  case  of  the  Bernstein  distribution: 


FT(t)  m  $ 


t-[D-  0]//i 
at/ p 


(2.8) 


Lu  and  Meeker  [20]  present  another  form  of  the  Bernstein  distribution  when  both  the 
intercept  and  the  slope  in  the  simple  linear  path  model  are  assumed  to  be  random, 
normally  distributed,  and  independent  of  each  other  (i.e.  0  is  replaced  by  01;  and 
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0  is  replaced  by  02  in  (2.5).)  In  this  case, 


FT(t)  ft  <f> 


t-\D  -  /ti]//x2  \ 


Another  nonlinear  model  presented  in  [20]  is 


rj(t)  =  0i  +  0exp(02f),02  >  0. 


(2.9) 


(2.10) 


Here,  0i  and  02  are  fixed,  and  the  resulting  time-to-failure  distribution  is  given  as 


FT(t)  ft 


t  -  [log (£>  -  0i )  -n\/<t> 2 
tr/02 


(2.n) 


In  this  case, 


T  ~  AT 


(\og(D  —  0i)  —  a* 

v  02  ’  0iy ' 


(2.12) 


Next,  Lu  and  Meeker  [20]  present  a  multivariate  normal  model  where  the  vector  of 
random  effects  0  follows  a  multivariate  normal  distribution  with  mean  vector  fig  and 
variance-covariance  matrix  Eg.  With  this  model,  the  time-to-failure  distribution  is 
expressed  as 


Fr(t)  —  FT(t',(f),  fig,  Eg).  (2-13) 

It  is  noted  in  [20]  that  there  is  generally  no  closed-form  expression  for  this  function. 
A  two-stage  method  for  estimating  the  random-effects  parameters  is  proposed  to 
avoid  the  computational  intensity  of  full  maximum  likelihood  estimation: 

1.  Fit  the  degradation  model  to  the  sample  path  for  each  sampled  unit  and  obtain 
the  Stage  1  estimates  of  the  model  parameters. 
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2.  Model  the  random-effect  parameters  with  a  (multivariate)  normal  distribution 
by  transforming  (if  necessary)  the  Stage  1  estimates. 

3.  Combine  the  transformed  Stage  1  estimates  of  the  model  parameters  to  produce 
estimates  of  </>,  pg,  and  Eg. 

Once  the  estimated  parameters  are  obtained,  Lu  and  Meeker  [20]  propose  a  method 
by  which  Monte  Carlo  simulation  is  used  to  generate  a  sufficiently  large  number  of 
random  sample  paths  from  the  assumed  path  model  (with  the  estimated  parameters). 
The  proportion  failing  is  then  used  as  a  function  of  time  as  an  estimate  of  Tr(f).  This 
method  is  particularly  useful  when  there  is  no  closed-form  expression  for  FT(t)  and 
when  numerical  transformation  methods  are  too  complicated.  The  authors  present 
the  following  algorithm: 

1.  Estimate  the  path-model  parameters  0,  /ig,  and  Eg; 

2.  Generate  N  simulated  realizations  0  of  0  from  iV(/i©,  E©)  and  obtain  the 
corresponding  N  simulated  realizations  0  of  0; 

3.  Compute  the  corresponding  N  simulated  failure  times  t  by  substituting  0  into 
T  =  t(0;  4>,  D,  77); 

4.  Estimate  Fxit)  =  (number  of  t  <  t)/N  for  any  desired  values  of  t. 

Lu  and  Meeker  [20]  then  apply  a  parametric  bootstrap  method  to  obtain  pointwise 
confidence  intervals  for  FT{t).  The  same  method  is  suggested  for  the  construction  of 
simultaneous  confidence  bands. 

Chao  [8]  examines  both  fixed  effect  and  random  effect  models.  A  simple  fixed 
effect  model  for  predicting  the  shelf  life  of  drugs  is  given  by 

F'ij  Oi{  T  fctij  T  ,  (2.14) 

where  YXJ  is  the  result  of  the  j-th  assay  of  the  i-th  batch  of  drugs,  t ij  is  of  dimension 
p  x  1  and  denotes  the  appropriate  regressor  (time),  and  (it  are  fixed  but  unknown 
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parameters,  and  etj  denotes  the  error  term.  In  this  case,  “fixed  effect  model”  refers 
to  the  random  selection  of  a  single  drug.  Chao  [8]  notes  that  cq  and  fa  can  be 
modelled  as  random  variables  for  random  effect  models,  which  refer  to  those  models 
considering  selection  of  a  batch  of  drugs  at  random  from  many  batches  at  hand. 
Chao  [8]  suggests  that  growth  curve  models  provide  a  more  general  framework  than 
the  determination  of  shelf  life.  In  cases  where  growth  is  bounded  above,  an  S-shaped 
curve,  or  “sigmoid”,  may  be  used  to  model,  from  beginning  to  end,  the  changing 
growth  rate.  Chao  [8]  proposes  the  following  model  for  growth: 


Vi  = 


a 


1  +  expfff  —  Bti 


+  hd  —  h  2, 


,  n. 


(2.15) 


This  model  is  often  called  a  logistic  growth  model. 

Chao  [8]  suggests  there  are  two  basic  approaches  used  to  describe  a  degradation 
process:  the  direct  approach  and  the  indirect  approach.  The  direct  approach  begins 
with  equations  and  uses  data  fitting  for  analysis.  The  indirect  approach  begins  with 
theoretical  justification  and  attempts  to  derive  more  theory  without  support  from 
data.  If  X  (t)  is  a  cumulative  damage  process  with 


X{t)  =  a  +  bt  +  W  (£), 


(2.16) 


where  W (t)  is  a  Brownian  motion,  then  T  has  a  distribution  which  is  inverse  Gaus¬ 
sian,  and  the  exact  solution  can  be  found  (cf.  Bhattacharyya  and  Fries  [5]). 

Chinnam  [9]  developed  a  degradation-based  approach  for  on-line  drill-bit  reli¬ 
ability  determination  through  the  use  of  neural  networks  for  modeling  degradation 
measures  and  self-organizing  maps  for  modelling  degradation  variation.  In  [9]  Chin¬ 
nam  proposes  a  failure  time  distribution  based  upon  a  nonlinear  regressive  degra¬ 
dation  model  with  predictions  of  future  degradation  levels  provided  by  the  artificial 
neural  network. 
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Crk  [11]  estimated  degradation  model  parameters  for  every  single  unit,  then 
combined  them  to  obtain  population  parameters  at  every  stress  level.  Crk  [11]  then 
uses  multiple  multivariate  regression  to  extrapolate  parameters  at  use  stress  levels. 
System  reliability  can  be  estimated  using  large  numbers  of  generated  failure  times 
derived  by  generating  large  numbers  of  model  parameters. 

Wu  and  Tsai  [28]  considered  the  case  where  a  few  of  the  degradation  paths  in 
a  test  differ  from  the  rest.  After  modelling  the  degradation  paths  via  nonlinear  re¬ 
gression,  the  authors  used  an  optimal  fuzzy  clustering  method  to  improve  estimation 
of  the  random  parameters  and  failure-time  distribution. 

Yacout,  et  al.  [29]  studied  a  mixed-effects  degradation  model  for  the  case 
of  cladding  strain  in  irradiated  fuel  pins.  Yacout,  et.al.  [29]  used  Monte  Carlo 
simulation  to  obtain  point  estimates  for  the  distribution  and  its  confidence  intervals 
of  this  mixed-effects  model  which  is  based  on  the  work  of  Lu  and  Meeker  [20] . 

Lawless  [18]  gave  an  overview  of  degradation  models  and  failure-time  distribu¬ 
tions,  including  shared  random  effects  models.  The  author’s  overview  also  includes 
Wiener  and  gamma  processes  and  references  the  random  growth  curve  in  [20]. 

As  this  brief  review  of  the  literature  attempts  to  show,  many  probabilistic  and 
statistical  methods  for  degradation-based  reliability  assessment  are  available  and 
can  be  extremely  useful  to  the  practitioner.  Although  many  statistical,  degradation- 
based  reliability  techniques  may  provide  accurate  lifetime  estimations,  the  use  of 
techniques  such  as  nonlinear  regression  requires  observation  and  measurement  of 
the  system’s  degradation  over  time  in  order  to  construct  a  failure  time  model  and 
estimate  its  parameters.  Frequently,  such  measurements  are  prohibitively  expensive 
or  even  impossible  to  obtain. 

Probabilistic  models  in  which  the  wear-inducing  environment  is  stochastically 
modelled  can  also  prove  to  be  useful  tools  for  reliability  assessment.  Unlike  statisti¬ 
cal,  degradation-based  models,  a  probabilistic  approach  does  not  require  degradation 
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measurements  during  the  lifetime  of  the  item.  However,  restrictions  on  either  the 
random  environment  or  wear  processes,  such  as  the  requirement  for  exponentially 
distributed  state  holding  times  seen  in  [16],  limit  the  scope  of  reliability  problems  to 
be  studied. 

In  this  thesis,  a  probabilistic  approach  discussed  by  Singpurwalla  [25]  and  thor¬ 
oughly  detailed  by  Kharoufeh  [16]  is  chosen  to  evaluate  the  failure  time  distribution 
of  a  system  under  the  influence  of  a  dynamic  environment.  The  objective  is  to  de¬ 
velop  a  lifetime  estimation  model  that  does  not  rely  upon  degradation  measurements 
during  the  system’s  lifetime  and  makes  only  minimal  assumptions  about  the  system’s 
operating  environment.  Such  a  model  provides  an  extremely  useful  tool  for  the  relia¬ 
bility  assessment  of  systems  operating  and  degrading  in  dynamic  environments.  The 
next  chapter  utilizes  and  relaxes  the  assumptions  of  Kharoufeh  [16]  and  Kharoufeh 
and  Sipe  [17]  and  provides  the  main  results  of  this  research  effort. 
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3.  Formal  Mathematical  Model 

This  chapter  provides  the  mathematical  model  used  to  assess  the  failure  time 
distribution  and  moments  of  failure  time  of  a  single-unit  system  experiencing  nonde¬ 
creasing  wear  due  to  a  multi-state  random  environment  process.  Since  it  is  assumed 
that  the  evolution  of  the  environment  depends  only  upon  its  current  state,  a  discrete- 
time  Markov  chain  (DTMC)  underlies  the  environment  process.  Furthermore,  since 
the  distribution  of  random  sojourn  times  is  assumed  to  be  nonexponential  or  even 
unknown,  the  environment  must  be  characterized  as  a  semi-Markov  process  (SMP). 
The  model  proposed  herein  converts  the  SMP  environment  process  to  a  CTMC  envi¬ 
ronment  process,  allowing  for  the  results  of  Kharoufeh  [16]  and  Kharoufch  and  Sipe 
[17]  to  be  applied. 

In  order  for  the  environment  to  be  characterized  as  a  SMP,  its  true  state  space 
must  first  be  partitioned  into  K  distinct  states.  Such  a  partitioning  may  be  achieved 
by  observation  and  measurement  or  by  consultation  with  subject  matter  experts. 
Having  partitioned  the  environment,  the  rate  at  which  each  state  of  the  environment 
causes  degradation  to  the  system  must  be  obtained.  In  real-world  applications,  it 
is  likely  that  subject  matter  experts  can  provide  insight  into  the  degradation  rates 
associated  with  the  states  of  a  particular  environment  on  a  chosen  system.  In  the 
proposed  model,  the  degradation  rate  associated  with  each  environment  state  is 
assumed  to  be  a  constant.  Next,  the  environment  process  is  observed  over  a  long 
period  of  time.  Holding  times  in  each  of  the  environment  states,  as  well  as  the 
numbers  of  transitions  between  states,  are  recorded.  Finally,  the  model  is  applied  to 
obtain  a  reliability  estimate  of  the  system  based  on  its  failure  time  distribution  and 
moments  of  failure  time.  A  detailed  discussion  of  this  process  follows  in  section  3.4. 

Due  to  their  relevance  to  this  thesis  research,  the  main  results  of  [16]  are 
reviewed  here.  Kharoufeh  [16]  derived  the  failure  time  distribution  and  moments  of 
a  single- unit  system  whose  cumulative  damage  over  time  is  a  continuous  wear  process, 
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{X(t)  :  t  >  0},  dependent  on  an  external  environment  process.  The  environment 
process  was  assumed  to  be  a  temporally  homogeneous,  Markov  process,  { Z{t )  :  t  > 
0}  on  a  a  finite  state  space  S  =  {1  Kharoufeh  [16]  derived  the  cumulative 

distribution  function  of  Tx,  the  time  until  failure,  by  using  transform  methods  to 
solve  a  first-order  linear  partial  differential  equation  satisfied  by  the  joint  probability 
distribution  of  the  Markov  additive  process  {(X(f),  Z(t))  :  t  >  0},  conditioned  upon 
the  initial  state  of  the  environment.  Using  the  derived  failure  time  distribution, 
Kharoufeh  [16]  presented  the  Laplace-Stieltjes  transform  of  the  failure  time  moments. 
Define, 


G(x,t )  :=  P{TX  <  t} 

as  the  unconditional  distribution  of  TXl  and  the  K  x  K  matrix 

V(x,t)  =  [Vi  ,j(x,t)], 


(3.1) 


(3.2) 


where 


Vi,j  (x,t)  =  P{X(t)  <  x,  Z(t)  =  j  \  Z( 0)  =  i}  (3.3) 

is  the  joint  probability  that,  at  time  t,  the  degradation  of  the  system  has  not  exceeded 
a  value  x,  and  the  environment  process  is  in  state  j  E  S  given  that  the  environment 
was  initially  in  state  j  G  5.  When  r  :  S  M+  is  defined  as  a  nonnegative  function, 
the  transform  of  the  joint  probability  above  is  given  as 

V  (n,  s)  =  [uKd  +  si  —  Q]-1,  (3.4) 

where  Re[u)  >  0,f?e(s)  >  0,  Rd  =  diag(r(l), ...,  r(K)),  and  Q  =  [qtJ]  is  the  in¬ 
finitesimal  generator  matrix  of  the  Markov  process,  {Z(t)  :  t  >  0}.  The  failure  time 
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distribution  of  a  single-unit  system  in  a  Markovian  environment  is  given  by 

G*(u,  s)  = - aV  (u,  s)l,  Re(u )  >  0,  Re(s)  >  0,  (3.5) 

s 

where  a  :=  [cq]  is  the  initial  distribution  vector  with  cq  :=  P{Z{ 0)  =  i},  and  1 
is  a  X- dimensional  column  vector  of  ones.  Kharoufeh  and  Sipe  [17]  simplified  this 
expression  to  a  one-dimensional  Laplace-Stieltjes  transform  with  respect  to  t,  and 
their  main  result  for  the  unit’s  failure  time  is  given  by 

Gx(s)  =  aexp(R^1(Q  —  sl)x)l.  (3.6) 

Let  £n(x)  :=  E[T™]  be  the  unconditional  nth  moment  of  the  failure  time.  Then, 
if  Z  has  initial  probability  distribution  a,  the  Laplace-Stieltjes  transform  of  the 
unconditional  nth  moment  of  the  failure  time  is  given  by 

=  n!a(tiRD  -  Q)“"l.  (3.7) 

Any  one  of  several  one-dimensional  numerical  Laplace  transform  inversion  algorithms 
can  be  used  to  invert  this  transform.  In  [16]  and  [17],  the  authors  provide  a  simplified 
method  for  estimating  the  reliability  of  a  system  based  solely  on  the  degradation 
effects  of  the  operating  environment  on  the  system,  where  the  environment  process 
is  assumed  to  be  characterized  as  a  CTMC. 

The  reliability  model  of  Equation  (3.6)  strictly  assumes  that  the  holding  times 
in  each  environment  state  are  exponentially  distributed  random  variables.  This 
requirement  facilitates  the  analysis  of  the  wear  process  via  a  CTMC  environment 
process.  However,  the  real  operating  environment  of  a  unit  may  not  necessarily  pos¬ 
sess  exponentially  distributed  state  holding  times.  This  thesis  proposes  a  method 
that  permits  the  relaxation  of  the  requirement  for  exponentially  distributed  holding 
times  while  still  applying  the  distribution  result  of  Equation  (3.6).  In  fact,  the  pro- 
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posed  technique  does  not  require  knowledge  of  the  distributions  of  the  state  holding 
times-only  real  observations  of  the  sojourn  times  in  each  state.  Furthermore,  since 
arbitrary  distributions  can  be  represented  exactly  or  approximately  by  phase-type 
distributions,  such  approximations  will  serve  as  the  key  to  the  relaxation  of  the  re¬ 
quirement  for  exponentially  distributed  state  holding  times  [23].  By  approximating 
nonexponential  state  holding  times  with  phase-type  distributions,  the  Markovian 
property  may  be  retained,  thereby  permitting  analysis  via  CTMC-based  techniques. 
In  the  following  subsections,  the  rudimentary  concepts  of  phase-type  distributions 
are  reviewed,  including  a  few  of  their  applications  and  techniques  for  obtaining  their 
representations. 

3.1  Phase-Type  Distributions 

In  general,  a  phase-type  distribution  can  be  described  as  some  mixture  of  ex¬ 
ponential  distributions,  each  representing  a  phase,  with  or  without  the  same  rate 
parameter.  The  phase-type  random  variable,  described  by  the  total  time  spent 
traversing  the  exponential  phases,  is  equivalent  to  the  time  to  absorption  of  an  un¬ 
derlying  Markov  process.  For  example,  a  fc-Erlang  random  variable,  the  sum  of  k 
independent  and  identically  distributed  exponential  random  variables,  is  equivalent 
to  the  absorption  time  of  an  underlying  /c-state  Markov  process.  As  noted  by  Perros 
[23],  the  value  of  phase-type  distributions  lies  both  in  their  possessing  the  Marko¬ 
vian  (memoryless)  property  and  their  usefulness  in  either  exactly  or  approximately 
representing  arbitrary  distributions.  These  two  characteristics  form  the  basis  for  the 
use  of  phase-type  distributions  in  this  work. 

The  use  of  phase-type  distributions  to  model  nonexponential  holding  times  is 
thoroughly  studied  in  the  literature.  One  useful  application  of  phase-type  distribu¬ 
tions  is  in  approximating  channel  holding  time  distributions  in  mobile  communica¬ 
tions  networks.  Jayasuriya,  et  al  [15]  used  generalized  Erlang  distributions  to  model 
channel  holding  time  distributions,  while  Ro  and  Trivedi  [24]  used  three-phase  hy- 
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poexponential  distributions.  Boucherie  and  Van  Dijk  [6]  apply  mixtures  of  Erlang 
distributions  as  approximations  to  call  length  distributions  and  the  distribution  of 
the  time  spent  in  a  cell.  Fang  and  Chlamtac  [13]  utilized  Erlang  and  hyper-Erlang 
call  holding  time  distributions  in  their  analysis  of  handoff  probability.  Phase-type 
distributions  have  also  been  used  to  approximate  service  times  at  different  stations 
in  a  production  line.  Viclalis  and  Papadopoulos  [26]  studied  the  transition  matrices 
of  production  lines  by  assuming  a  two-phase  Coxian  distribution  for  service  times, 
where  the  Erst  phase  describes  the  service  period,  and  the  second  phase  describes 
the  repair  period  which  is  reached  with  some  specified  probability. 

The  following  is  a  slight  modification  of  the  description  of  the  representation  of 
a  phase- type  distribution  as  found  in  [23] .  Recall  that  the  total  time  spent  traversing 
the  k  phases  of  a  phase-type  distribution  is  equivalent  to  the  time  to  absorption  of 
some  underlying  Markov  process.  Any  phase-type  distribution  can  be  represented  by 
a  triplet  (/3,  T.  T°),  where  [3  =  (/3,  0)  is  the  1  x  (k  +  1)  initial  probability  vector  of  the 
overall  phase  space  of  the  underlying  Markov  process,  (3  is  the  1  xk  vector  giving  the 
probabilities  that  the  underlying  Markov  process  will  start  in  phase  i,i  =  1,  2, ...,  k, 
T  is  the  k  x  k  matrix  of  transition  rates  among  the  Erst  k  phases,  and  T°  is  the  kx  1 
vector  of  transition  rates  out  of  the  first  k  phases  into  the  absorbing  phase  k  +  1. 
Note  that  (3  =  (/3,  0)  implies  that  the  underlying  Markov  process  never  begins  in 
the  absorbing  phase  k  +  1.  The  rate  matrix  of  the  phase-type  distribution  (and  the 
underlying  Markov  process)  is 


S 


(3.8) 


Using  the  representation  (/3,  T.  T°),  the  density  function  of  a  phase-type  distribution 
is  given  by 


f(x )  =  (3 ex p  (Tx)T°,x  >  0, 


(3.9) 
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with  Laplace  transform 


f*{s)  =  (3{sl  -  T)_1T°,  Re(s)  >  0.  (3.10) 

The  nth  moment  (n  >  1)  of  a  phase-type  distribution  with  representation  (/3,  T)  is 
given  by 


E[Xn ]  =  (-l)nn\pT~nl.  (3.11) 

The  next  two  sections  describe  various  phase-type  distributions,  their  structures  of 
their  representations,  and  techniques  for  obtaining  the  parameters  of  their  represen¬ 
tations. 

3.2  Types  of  Phase- Type  Distributions 

In  this  subsection,  various  types  of  phase-type  approximations  are  briefly  ex¬ 
amined.  In  the  next  section,  the  technique  used  herein  to  construct  phase-type 
approximations  when  state  holding  times  are  known,  either  via  observed  data  or 
known  parametric  probability  distributions,  is  reviewed. 

3.2.1  General  Phase-Type  Distribution 

The  general  case  of  a  phase-type  distribution,  denoted  by  PHkl  is  one  in  which 
the  exponential  phases  are  not  necessarily  visited  sequentially.  In  other  words,  a 
visit  to  any  phase  i,i  —  1  can  be  followed  by  a  visit  to  any  other  phase 

j,j  =  1  ,...k,i  7^  j,  with  probability  prior  to  absorption.  It  is  assumed  here 
that  absorption  occurs  upon  reaching  an  unseen  absorbing  phase.  Furthermore, 
each  phase  i,i  —  1,2 has  exponential  rate  parameter  .  Figure  3.1  gives  a 
graphical  depiction  of  the  general  case  of  a  phase-type  distribution.  For  the  sake  of 
clarity,  some  of  the  transition  probabilities  have  been  omitted. 
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Figure  3.1  Graphical  depiction  of  a  /c-phase,  phase-type  distribution. 

3.2.2  Coxian  Distribution 

The  k- phase  Coxian  distribution,  denoted  by  C*,  consists  of  a  sequence  of  k 
exponential  distributions  with  respective  rate  parameters  /q ,  i  —  1,  2, ...,  k.  Asa  spe¬ 
cial  case  of  the  general  phase-type  distribution,  the  Coxian  distribution  is  sequential 
in  its  phases  and  does  not  permit  transitions  from  phase  j  to  phase  i  when  j  >  i. 


► 


Figure  3.2  Graphical  depiction  of  a  fc-phase  Coxian  distribution. 

In  Figure  3.2,  <q,  i  —  0, 1, ...,  k  *-  1,  denotes  the  probability  of  departing  phase  i  and 
entering  phase  i  +  1,  while  bi,  i  =  0, 1, ...,  k  —  1,  denotes  the  probability  of  departing 
phase  i  and  entering  the  absorbing  phase.  Also,  note  that  ao  denotes  the  probability 
of  entering  phase  1,  while  b0  denotes  the  probability  of  zero  service  time.  A  useful 
property  of  the  Coxian  distribution  is  that  it  can  exactly  represent  any  distribution 
having  a  rational  Laplace  transform  [23].  Moreover,  Perros  [23]  gives  the  Laplace 
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transform  of  a  Coxian  distribution  as 


k  i 

P(s)  =  bo  +  V'ao-.  .a, _i6, TT  — ^ —  (3.12) 

£r  s + ft 

where  a*  +  6,:  =  1,  f  =  1,  2, k  —  1,  and  =  1. 

3.2.3  Erlang  Distribution 

The  Erlang  distribution  is  a  special  case  of  the  Coxian  in  which  each  exponen¬ 
tial  phase  has  equal  rate  parameter  /q  and  there  is  zero  probability  of  entering  the 
absorbing  phase  (or  any  phase  other  than  phase  i  +  1)  until  the  kth  phase  is  reached 
[23].  Figure  3.3  gives  a  graphical  depiction  of  this  distribution. 


Figure  3.3  Graphical  depiction  of  a  k- phase  Erlang  distribution  with  rate  parameter  g. 

A  more  general  form  of  the  Erlang  distribution,  depicted  in  Figure  3.4,  is  obtained 
by  allowing  distinct  exponential  phase  rate  parameters  q*,  i  —  1,  2, ...,  k. 


Figure  3.4  Graphical  depiction  of  a  /c-phase  Erlang  distribution  with  distinct  rate  pa¬ 
rameters. 

The  generalized  Erlang  distribution  is  a  /c-phase  Erlang  distribution  that  transitions 
from  phase  1  to  phase  2  with  probability  a  and  from  phase  1  to  the  absorbing  phase 
with  probability  1-a.  If  phase  2  is  reached,  then  all  phases  after  phase  2  are  also 
reached  (sequentially).  Note  that  an  alternative  definition  of  the  generalized  Erlang 
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distribution,  equivalent  to  the  multiple  rate  Erlang  distribution  described  above,  is 
implied  in  [21]. 


Figure  3.5  Graphical  depiction  of  a  fc-phase  generalized  Erlang  distribution  with  single 
rate  parameter. 


Since  the  time  required  to  traverse  all  phases  of  a  /e-phase  Erlang  distribution  is 
the  sum  of  k  independent  and  identically  distributed,  exponential  random  variables 
with  common  rate  parameter  /i,  the  Laplace  transform  (LT)  of  the  A;- phase  Erlang 
distribution  is  given  by 


rw  = 


(3.13) 


and  the  LT  of  the  fc-phase  Erlang  distribution  with  distinct  rate  parameters  /.q ,  i  = 
1,  2,  •••,  k,  is 


/*(*) 


(3.14) 


3.3  Phase-type  Approximation  Technique 

In  this  section,  the  procedure  used  in  this  thesis  to  approximate  observed  data 
or  parametric  probability  distributions  with  phase-type  distributions  is  presented. 
Table  3.1  provides  some  guidelines  for  choosing  the  particular  phase- type  approxi¬ 
mation  to  use  based  on  the  squared  coefficient  of  variation  of  the  observed  data  or 
parametric  probability  distribution  [23].  More  specifically,  suppose  X  is  an  arbitrary 
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nonnegative  random  variable.  The  objective  is  to  obtain  a  phase-type  distribution 
to  approximate  the  true  distribution  of  X.  The  squared  coefficient  of  variation  for 
the  random  variable  X  is 


2  Var  [X] 

cdvF 


(3.15) 


Table  3.1  Guidelines  for  selection  of  phase- type  approximation. 
Range  of  c2  Distribution 

c2  >  1  three-moment,  2-phase  Coxian 

0.5  <  c2  <  1  two-moment,  2-phase  Coxian 

c2  <  0.5  two- moment,  /e-phase  generalized  Erlang 


The  following  calculations  follow  from  [4]  and  [23]. 


3.3.1  Approximation  via  three-moment,  2-phase  Coxian  distribution 

When  c2  >  1,  the  observed  holding  time  data  or  parametric  probability  dis¬ 
tribution  is  approximated  with  a  three-moment,  2-phase  Coxian  distribution.  The 
Laplace  transform  (LT)  of  this  approximation  is  given  as 

y*/  x  _  /xis(l  —  a)  +  P1P2 

S2  +  (/ii  +  fl2)s  +  /Ui/i2 

where  s  is  a  complex  transform  variable,  and  /Ui,/i2,  and  a  denote  the  parameters 
of  the  Coxian  representation  as  described  in  section  3.2.  After  taking  successive 
derivatives  of  the  above  LT  and  evaluating  at  s  —  0,  the  first  three  moments  of  the 
C'2  distribution  may  be  obtained.  Let  Y  be  a  2-phase  Coxian  distributed  random 
variable.  Estimates  (or  exact  values)  of  the  first  three  moments  (mi,  m2,  and  m3, 
respectively)  of  the  observed  data  or  parametric  probability  distribution  are  set  equal 
to  the  three  C'2  moments  and  are  expressed  as  (cf.  [4]) 


.  ,  1  a 

mi  =  E[Y]  = - 1 - , 

H 1  H  2 


(3.17) 
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(3.18) 


m2  =  E[Y2]  = 

H'l 


2  2(1  -  a)  2ajiiH2  -  2 a(ni  +  //2)2 


f4t4 


and 


m,  =  E[Y3)  =  3 


3  _  6(1  -  a)  12a/ii/i2(/ii  +  /i2)  -  6a(/ii  +  /i2)3 


(3.19) 


For  convenience,  the  following  equalities  are  introduced: 


A  —  Hi  +  /i2 


(3.20) 


and 


i?  —  /ii/i2. 


(3.21) 


Upon  substitution  and  simplification  of  Equations  (3.17),  (3.18),  and  (3.19),  the 
following  expressions  are  obtained: 


.  1  m2B 

A  = - 1 - 

m  i  2m  i 


(3.22) 


and 


B  = 


12ml  ~  6'm2 

3/77.2  —  2/771/773 


(3.23) 


Following  [4],  the  parameters  of  the  three-moment,  2-phase  Coxian  distribution  are 


// 1  —  A  + 


Vdl2  -  AB 


(3.24) 


fi2  —  A  —  fj>i, 


(3.25) 
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and 


(3.26) 


a  =  —(mini  —  1). 

Pi 

The  three-moment,  2-phase  Coxian  distribution  has  the  representation 


T  = 


(3.27) 


>9  =  (1,0,0), 


(3.28) 


T°  —  [(1  —  a)p  i ,  /x2] . 


(3.29) 


3.3.2  Approximation  via  two-moment,  2-phase  Coxian  distribution 


When  0.5  <  c2  <  1,  a  two-moment,  2-phase  Coxian  dsitribution  is  used  to 
approximate  the  observed  holding  time  data  or  parametric  probability  distribution. 
Following  a  moment-matching  algorithm  similar  to  that  for  the  three-moment  2- 
phase  Coxian  distribution,  the  parameters  of  the  two-moment  variant  of  the  2-phase 
Coxian  are 

pi  =  — ,  (3.30) 

m  i 

P'2  —  2’  (3.31) 

m\Cz 


and 


(3.32) 


The  two-moment,  2-phase  Coxian  distribution  has  the  same  representation  as  the 
three-moment,  2-phase  Coxian  distribution  described  above.  It  is  important  to  note 
that  using  either  variant  of  the  2-phase  Coxian  distribution  as  an  approximation  to 
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observed  holding  time  data  or  a  parametric  probability  distribution  helps  to  minimize 
the  growth  of  the  overall  state  space.  The  drawback  of  state  space  growth  is  further 
discussed  in  section  3.4.1. 


3.3.3  Approximation  via  two-7noment,  k-phase  Erlang  distribution 


When  the  squared  coefficient  of  variation  is  less  than  0.5,  the  observed  holding 
time  data  or  parametric  probability  distribution  is  approximated  by  a  two-moment, 
/e-phase  generalized  Erlang  approximation.  Following  [23],  the  number  of  phases  k 
is  chosen  such  that 


1 

k 


<  c2  < 


1 

V~i~y 


(3.33) 


In  the  numerical  examples  of  chapter  4,  a  simple  linear  program  is  implemented  to 
approximate  the  minimum  number  of  phases  k.  The  remaining  parameters,  p  and 
a ,  corresponding  to  those  shown  in  Figure  3.5,  are  computed  as 


2  kc2  +  k-2-(k2 +  4-  4kc 2)1/2 
2(c2 +  !)(&  —  1) 


(3.34) 


and 


1  +  [k  —  l)a 
m\ 


The  /c-phase  generalized  Erlang  distribution  is  then  represented  by 


T 


^  —  p  ap  0  0 

0  —  p  p  0 

0  0  —  p  p 

0  0  0  —  p 


0  N 
0 

0 


y  o  o  o 


••  v  \ 

o  -p  ) 


/3=(  1,0,  ••  •  ,0) 


(3.35) 


(3.36) 


(3.37) 
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T°  =  [(l-a)/^,0, ...  ,0,/i]. 


(3.38) 


3-4  Phase-Type  Approximations  for  Degradation- Based  Reliability 

In  real-world  situations,  when  analyzing  the  wear  process  of  a  single-unit  sys¬ 
tem  subject  to  a  multi-state  stochastic  environment  process,  knowledge  of  the  en¬ 
vironment  process  will  likely  be  limited  to  the  degradation  rate  associated  with 
each  state  and  observations  of  the  state  holding  times.  Since  there  is  no  guarantee 
that  these  holding  times  will  be  exponentially  distributed,  the  environment  process 
cannot  be  assumed  to  be  a  continuous-time  Markov  chain.  However,  phase-type  dis¬ 
tributions,  such  as  those  described  in  section  3.2,  may  be  used  as  approximations  to 
these  nonexponential  distributions  to  impose  the  Markovian  property.  That  is,  the 
approximated  environment  process  may  be  modelled  as  a  continuous-time  Markov 
chain  (CTMC),  and  subsequently  a  new  (and  expanded)  infinitesimal  generator  ma¬ 
trix  may  be  used  directly  in  Equation  (3.6). 

3.4-1  Conversion  of  the  Environment  Process  to  a  CTMC 

Consider  a  single-unit  system  subject  to  wear  over  time  due  to  a  multi-state 
semi-Markov  environment  process,  { Z{t )  :  t  >  0},  with  finite  state  space  S  = 
{1,2 , ,K}.  For  the  sake  of  brevity,  the  SMP  is  abbreviated  as  Z.  The  objective 
is  to  convert  Z  into  a  CTMC.  To  accomplish  this,  each  sojourn  time  distribution 
Hj,  i  e  S,  of  Z  is  approximated  by  a  phase- type  distribution  Hi.  The  state  space  S  of 
the  resulting  CTMC  consists  of  all  of  the  phases  in  the  K  phase-type  approximations 
Hi. 

Initially,  the  environment  process  must  be  observed  over  a  long  time  interval  of 
length  r.  The  holding  times  in  the  K  states  and  the  random  number  of  transitions 
from  one  state  to  another  are  observed  and  recorded.  The  observed  rate  of  transition 
from  state  i  to  state  j,  i,j  €  S,  over  r  is  computed  and  denoted  as  cpy  The  K  x  K 
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transition  rate  matrix  Q,  whose  elements  are  ql3 ,  is  then  formed.  After  computing 
the  required  moments  of  the  observed  holding  times  in  each  state,  the  K  holding  time 
distributions  (Ht,  i  =  1,2 are  approximated  by  K  phase-type  distributions 
(Hi}  i  =  1,  2, ...,  K)  by  applying  the  techniques  of  section  3.3.  Recall  that  either  a  2- 
phase  Coxian  or  a  fc-phase  generalized  Erlang  distribution  is  chosen  to  approximate 
each  holding  time  distribution  depending  on  the  estimated  squared  coefficient  of 
variation  of  each  state’s  observed  holding  times. 

A  potential  drawback  of  this  technique  is  the  possibility  of  a  state  space  ex¬ 
plosion.  For  example,  consider  a  five-state  semi-Markov  environment  process  with 
arbitrary  sojourn  time  distributions  Hi,  i  =  1,  2, ...,  5.  If  each  distribution  is  approx¬ 
imated  by  a  four-phase,  generalized  Erlang  distribution,  the  resulting  CTMC  will 
have  twenty  states,  represented  by  the  phases  used  to  approximate  each  distribution. 
In  an  effort  to  minimize  the  dimensionality  of  the  new  state  space,  and  therefore  ease 
the  computational  burden,  one  should  use  techniques  that  result  in  a  minimal  num¬ 
ber  of  phases  for  each  sojourn  time  distribution.  The  recent  work  of  Osogami  and 
Harchol-Balter  [22]  aims  at  minimizing  the  number  of  phases  while  maximizing  the 
accuracy  of  the  phase-type  approximation. 

Define  ki  as  the  number  of  phases  used  to  approximate  sojourn  time  distri¬ 
bution  Hi.  For  any  state  i  G  S,  the  underlying  Markov  process  has  ki  +  1  phases, 
including  ki  exponential  phases  plus  one  absorbing  phase.  In  order  to  complete  the 
conversion  of  Z  into  a  CTMC,  the  infinitesimal  generator  matrix  Sh,  whose  elements 
represent  the  transition  rates  among  all  phases  of  the  newly  partitioned  K  states, 
must  be  created.  Since  the  initial  probability  vector  /3j  of  each  phase-type  approx¬ 
imation  is  of  the  form  (1,  0, . . . ,  0),  transitions  from  state  i  to  state  j  for  i  ^  j  are 
limited  to  transitions  from  the  absorbing  phase  of  state  i  to  the  Erst  phase  of  state 
j.  In  other  words,  any  state  must  be  exited  via  its  absorbing  phase,  and  a  new  state 
must  be  entered  via  its  first  phase. 
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As  an  example,  suppose  one  seeks  to  convert  (approximately)  a  instate  semi- 
Markov  environment  process  into  a  CTMC  using  2-phase  Coxian  distributions  for 
each  of  the  K  states.  After  the  number  of  transitions  among  the  K  states  of  the 
original  environment  process  have  been  observed  and  recorded,  the  observed  rates  of 
transition  from  state  i,i  =  1,2, ,  K ,  to  state  j,j  =  1,  2, . . . ,  K,  are  computed  as 


.  Nr(i,j) 

9'3  HT(i ) 


(3.39) 


where  NT(i,j)  is  the  total  number  of  transitions  from  state  i  to  state  j  in  time 
r,  and  HT{i)  is  the  total  time  spent  in  state  i  during  time  r.  Figure  3.6  gives  a 
graphical  depiction  of  the  approximation  of  such  an  environment  when  K  =  3.  This 


State  1 


Figure  3.6  Graphical  depiction  of  a  3-state  environment  process  with  2-phase  Coxian 
sojourn  time  approximations. 

Jl -state  CTMC  environment  process  possesses  an  infinitesimal  generator  matrix  of 
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the  following  form: 


where 


and  Tj  and  T° 
tivcly,  of  state  i 


( 


^11 

*12  •• 

■  *1 K 

^21 

*22  •• 

■  *2  K 

^  K\ 

*  K2  ■  ■ 

■  ^I<K 

T 

rpo 

i 

\ 

\ 


o  <hM 


,  ieS, 


J 


( 


\ 


(3.40) 


(3.41) 


0  0  0 

*o  =1  0  0  0  ,  i^j  es,  (3.42) 

y  qijM  0  0  y 

are  the  phase-type  representative  T  matrix  and  T°  vector,  respec- 
(see  section  3.1). 


The  elements  of  i  e  S,  represent  the  rates  of  transition  among  the  kt  +  1 
phases  of  the  newly  partitioned  state  i ,  and  the  elements  of  i  ^  j  G  S,  represent 
the  rates  of  transition  from  the  kt  +  1  phases  of  state  i  to  the  kj  +  1  phases  of  state 
j.  In  other  words,  \E 'a,i  G  S,  gives  the  transition  rates  among  the  ki  +  1  phases  of 
a  single  state,  and  \E 'ij,i  ^  j  e  S,  gives  the  transition  rates  from  the  ki  +  1  phases 
of  one  state  to  the  kj  +  1  phases  of  another  state.  It  is  helpful  to  note  that,  in 
the  general  case  where  /e-phase,  phase-type  distributions  are  used  to  approximate  K 
environment  states,  the  dimensions  of  are  (ki  +  1)  x  (ki  + 1),  and  the  dimensions 
of  i  ^  j,  are  (kt  + 1)  x  (kj  +  1).  It  is  also  important  to  note  that  the  holding  time 
in  the  absorbing  phase  of  any  state  is  theoretically  instantaneous,  thereby  implying 
an  infinite  transition  rate.  Therefore,  in  numerical  implementations,  it  is  necessary 
that  the  transition  rate  from  the  absorbing  phase  of  any  state  to  either  itself  or  the 
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first  phase  of  any  other  state  is  multiplied  by  a  very  large  number  M.  In  other  words, 
each  <jij  in  th  is  multiplied  by  M  to  create  a  very  large  transition  rate. 

Given  that  the  wear  rate  in  environment  state  i  is  r(i),  i  =  1,2,... ,  K,  the 
wear  rate  matrix  of  the  semi-Markov  environment  process  is 


R 


D  ~ 


(  r(l)  0  ...  0 

0  r( 2)  i 

;  ■•.  ••.  o 

y  0  ...  0  r(K)  J 


(3.43) 


Having  expanded  the  original  generator  matrix  to  account  for  transitions  among 
the  phases  of  the  approximated  environment  states,  the  corresponding  wear  rate 
matrix  must  also  be  expanded.  The  wear  rate  r(i)  of  each  environment  state  i, 
i  =  1,2, .. . ,  K ,  is  assumed  to  be  the  same  for  all  kt  +  1  phases  of  the  phase-type 
approximation  of  state  i.  Therefore,  if  each  state  is  approximated  via  a  2-phase 
Coxian  distribution,  the  expanded  wear  rate  matrix  of  the  K -state  example  is  given 
as 


(  r(l)  0  0  0  0 

0  r(l)  000 
0  0  r(l)  0  0 

0  0  0  r(2)  0 

0  0  0  0  r(  2) 

0  0  0  0  0 

0  0  0  0  0 

0  0  0  0  0 

v  0  0  0  0  0 


0  ...  0  0  0  ^ 

0  ...  0  0  0 

0  ...  0  0  0 

0  ...  0  0  0 

0  ...  0  0  0 

r(2) 

'-.0  0  0 
...  0  r(K)  0  0 

...  0  0  r(K)  0 

...  0  0  0  r(K)  j 


(3.44) 
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In  the  next  subsection,  th  and  A  are  used  to  compute  the  failure  time  distribution 
and  moments  of  failure  of  the  single-unit  system  under  study. 

3-4-2  Computation  of  Failure  Time  Distribution  and  Moments 

Once  the  semi-Markov  environment  has  been  converted  into  a  CTMC  with 
infinitesimal  generator  matrix  \I/  and  wear  rate  matrix  A,  the  results  derived  by 
Kharoufeh  [16]  and  Kharoufeh  and  Sipe  [17]  may  be  applied  to  compute  the  cumu¬ 
lative  distribution  function  of  the  failure  time  of  the  single-unit  system  under  study. 
Figure  3.7  provides  a  summary  of  the  overall  process.  Applying  Equations  (3.6)  and 
(3.7),  the  Laplace-Stieltjes  transform  of  the  failure  time  distribution  is  given  by 

Gx(s )  =  a;  exp(A~1(\Er  —  sl)x)l,  (3.45) 

and  the  Laplace-Stieltjes  transform  of  the  nth  moment  of  unit  failure  time  is  given 
by 

CM  =  n\a{uK  -  ^)~nl,  (3.46) 

where  a  is  the  1  x  K  initial  probability  vector  determining  the  beginning  state  of 
the  environment  process.  Any  one  of  several  one-dimensional  inversion  algorithms 
may  be  used  to  numerically  invert  these  transforms  to  obtain  Gx(t)  and  fn(x),  the 
failure  time  distribution  and  nth  moment  of  unit  failure  time,  respectively.  The 
next  chapter  provides  three  numerical  examples  that  illustrate  the  conversion  of  a 
semi-Markov  environment  process  to  a  CTMC  environment  process  via  phase-type 
approximations.  The  failure  time  distribution  and  lower  moments  of  failure  time  are 
estimated  for  each  system  under  consideration. 
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Figure  3.7 


Flow  diagram  for  reliability  assessment  in  a  SMP  environment. 
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4.  Numerical  Results 


In  this  chapter,  applications  of  the  reliability  evaluation  process  described 
in  chapter  3  are  illustrated.  Using  phase-type  approximation  techniques,  the  fail¬ 
ure  time  distributions  for  single-unit  systems  subject  to  semi-Markov  environment 
processes  will  be  evaluated  via  the  results  previously  derived  for  systems  subject  to 
environment  processes  characterized  as  CTMCs.  Additionally,  the  first  and  second 
moments  of  failure  time  will  also  be  evaluated.  The  results  are  directly  applicable 
to  real-world  situations. 

4-1  Degradation  of  Brake  Pad  Material 

Consider  a  new  brake  pad  installed  on  an  automobile  equipped  with  disc  brakes. 
In  a  disc  brake  system,  a  caliper  holds  the  brake  pad  and  pushes  it  against  a  rotor  on 
the  automobile’s  wheel  when  the  brake  pedal  is  depressed.  In  this  example,  assume 
that  the  brake  pad  is  always  in  contact  with  the  rotor.  (In  modern  disc  brake  systems, 
the  brake  pad  sits  clear  of  the  rotor  when  the  brakes  are  not  applied.)  When  the 
brakes  are  applied,  the  caliper  pushes  the  brake  pad  against  the  rotor  with  a  force 
proportional  to  the  pressure  applied  to  the  brake  pedal.  The  friction  generated 
between  the  brake  pad  and  the  rotor  causes  the  automobile  to  slow  down.  Over 
time,  the  brake  pad  experiences  normal  wear  due  to  the  friction  between  the  pad 
and  the  rotor  and  eventually  reaches  a  failure  threshold  x.  Reaching  this  threshold 
may  not  necessarily  lead  to  system  failure  but  may  be  used  to  indicate  the  need  for 
preventive  maintenance  of  the  braking  system.  The  degradation  rate  of  the  brake 
pad  is  linear  and  completely  determined  by  the  system’s  random  environment,  which 
consists  of  three  states.  State  1,  which  is  described  by  the  normal  contact  of  the 
brake  pad  with  the  rotor,  causes  the  brake  pad  to  wear  with  rate  r(l).  States  2 
and  3,  achieved  by  depressing  the  brake  pedal  with  two  distinct  pressures,  cause  the 
brake  pad  to  wear  with  rates  r(2)  and  r(3),  respectively.  It  is  assumed  that  the 
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system’s  random  environment  is  characterized  by  a  semi-Markov  process  with  state 
space  S  =  {1,2,3}.  Let  {Xn  :  n  >  0}  be  a  DTMC  embedded  within  the  SMP  such 
that  Xn  is  the  state  of  the  environment  just  after  the  nth  transition.  The  DTMC 
has  transition  probability  matrix 


P 


^  0  a  1  —  a  ^ 

b  0  1  —  6 

^  c  1  —  c  0  j 


where  a  is  the  probability  of  transitioning  from  state  1  to  state  2,  b  is  the  probability 
of  transitioning  from  state  2  to  state  1,  and  c  is  the  probability  of  transitioning  from 
state  3  to  state  1. 

Let  X(t )  denote  the  cumulative  degradation  of  the  brake  pad  at  time  t,  let 
Re  be  the  diagonal  matrix  of  wear  rates  r(l),  r( 2),  and  r(3),  and  let  Tx  denote  the 
random  time  required  for  the  brake  pad’s  cumulative  degradation  to  reach  amount 
x.  The  cumulative  distribution  function  of  Tx  is  desired. 

In  order  to  apply  the  results  of  [16]  for  the  failure-time  distribution,  the  semi- 
Markov  environment  process  is  first  converted  into  a  CTMC  via  phase-type  ap¬ 
proximations  of  the  state  holding  time  distributions.  After  observing  the  random 
environment  over  time  r,  the  transition  rates  i,j  G  S,  are  computed  for  the 
transitions  among  the  environment’s  three  states.  This  is  accomplished  by  dividing 
the  number  of  times  the  system  transitions  from  state  i  to  state  j,  i  ^  j,  by  the 
total  amount  of  time  the  system  spends  in  state  i  for  i,j  e  S.  Additionally,  the 
observed  holding  times  in  each  state  of  the  environment  are  recorded  and  used  to 
construct  phase-type  approximations  to  the  holding  time  distributions.  The  result¬ 
ing  representations  of  the  phase-type  distributions  and  the  observed  state  transition 
rates  are  used  to  construct  T',  the  generator  matrix  of  the  resulting  CTMC.  LIsing 
the  generator  matrix,  Tf,  the  expanded  degradation  rate  matrix,  A,  and  the  initial 
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probability  distribution  of  the  environment  states,  a,  as  input,  Equation  (3.45)  is 
applied  to  calculate  the  solution  of  the  failure  time  distribution  in  the  transform. 

For  this  example,  the  wear  rates  are  r(l)  =  0.1,  r( 2)  =  1.0,  and  r(3)  =  2.0. 
In  order  to  conduct  a  simulation  of  the  environment  process,  the  state  holding  time 
distributions  are  assumed  to  be  known  and  are  given  in  Table  4.1.  Additionally, 

Table  4.1  State  holding  time  distributions  for  the  brake  pad  example. 


State 

Distribution 

1 

Weibull(2,4) 

2 

Weibull(3,5) 

3 

Weibull(4,6) 

the  transition  probability  matrix  values  are  a  =  0.7,6  =  0.6,  and  c  =  0.2,  and  the 
brake  pad  damage  threshold  is  set  at  x  =  30.0.  Furthermore,  it  is  assumed  that  the 
environment  process  always  begins  in  State  1.  Therefore,  the  required  input  matrices 
are 


and 


/ 


P  = 


V 

/ 


R/5  — 


V 


a  = 


0.0 

0.7 

0.3 

0.6 

0.0 

0.4 

0.2 

0.8 

0.0 

0.1 

0.0 

0.0 

0.0 

1.0 

0.0 

0.0 

0.0 

2.0 

(l  o  o). 


\ 

/ 

\ 

/ 


J 


After  observing  the  system  for  r  =  10000  time  units,  the  T,  matrices  and  T°  vec¬ 
tors  of  the  phase-type  representations  of  the  three  holding  time  distributions  were 
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estimated  as 


Ti  = 


T2  = 


(  . 

-8.826 

8.562 

0.000 

0.000  ^ 

(  0.2633 

\ 

0.000  - 

-8.826 

8.826 

0.000 

0 

= 

rpio 

1  1  — 

i 

0.000 

0.000 

-8.826 

8.826 

0 

V 

0.000 

0.000 

0.000 

-8.826  ) 

y  8.8256 

) 

-15.2250 

15.1090 

0 

0 

0 

0 

0 

0 

0 

-15.2250 

15.2250 

0 

0 

0 

0 

0 

0 

0 

-15.2250 

15.2250 

0 

0 

0 

0 

0 

0 

0 

-15.2250 

15.2250 

0 

0 

0 

0 

0 

0 

0 

-15.2250 

15.2250 

0 

0 

0 

0 

0 

0 

0 

-15.2250 

15.2250 

0 

0 

0 

0 

0 

0 

0 

-15.2250 

15.2250 

0 

0 

0 

0 

0 

0 

0 

-15.2250 

rpO  _ 

±2  ~ 


V 


0.1135 

0 

0 

0 

0 

0 

0 

15.2290 


\ 


J 


T3  =  [Sij]  for  i  =  1, . . . ,  13,  j  =  1, . . . ,  13, 


where 


and 


-22.4030, 

22.3580, 

22.4030, 

0 


i  =  j 

i  =  l,j  =  2 

otherwise 


T3  =  hi]  for  *  =  13, 
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where 


0.0450, 

i  =  1 

o, 

*  =  2,...,  12 

22.4030, 

i  =  13 

Furthermore,  the  observed  transition  rates  among  the  three  environment  states  are 
presented  as  the  elements  of  Q,  where 


Q 


(  -2.2569  1.5809  0.6760  ^ 

1.1501  -1.9162  0.7661 

v  0.3403  1.3862  -1.7265  ) 


Using  the  above  phase- type  representations  and  the  qij  elements  of  Q,  the  28  x  28 
generator  matrix  of  the  newly  formed  CTMC  is 


^12 

^13  \ 

^21 

^22 

CO 

CM 

* 

V 

* 

CO 

^32 

^33  / 

where 


where 


where 


T\  T? 

ii  —  I  i  > 

0  -22569.0 


^12  =  [Aj\  for  i  =  1, , . , ,  5,  j  =  1, . . . ,  9, 


Ip ij 


15809.0,  i  =  5,  j  =  1 


0  otherwise 
^13  ={fpij\  for  *  =  1,  -  -  - ,  5,  j  =  1,  —  ,  14, 


'P’ij  — 


6760.2,  i  =  5,  j  =  1 
0  otherwise 
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^21  =  [Aj]  for  i  =  1,  •  •  • ,  9,  j  =  1, . . . ,  5, 


where 


Aj  = 


11501.0,  i  =  9,  j  =  1 
0  otherwise 


where 


where 


22 


0  -19162.0 

/ 

^23  =  [Aj]  for  *  =  1,  •  •  •  $  9,  j  =  1, . . . ,  14, 


7660.9,  i  —  9,  j  —  1 

A3  =  { 

0  otherwise 

^3i  =  [Aj\  for  i  =  1,  •  •  • ,  14,  j  =  1, . . . ,  5, 


Aj  = 


3403.1,  i  =  U,j  =  l 


0  otherwise 

^32  =  [Aj]  for  i  =  1,  •  •  - ,  14,  j  =  1, . . . ,  9, 


where 


and 


Aj 


13862.0,  i  =  U,j  =  l 
0  otherwise 


^33  = 


t3  t° 

0  -17265.0 


As  noted  in  chapter  3,  the  holding  time  in  the  absorbing  phase  of  any  state 
is  theoretically  instantaneous,  thereby  implying  a  theoretically  infinite  transition 
rate.  Therefore,  the  transition  rate  from  the  absorbing  phase  of  any  state  to  either 
itself  or  the  first  phase  of  any  other  state  is  multiplied  by  a  large  number  M  when 
constructing  the  infinitesimal  generator  matrix  Tr.  In  this  case,  and  in  subsequent 
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examples,  M  =  10000.  The  28  x  28  expanded  degradation  rate  matrix  is 


A 


(  An  0  0  ^ 

0  A22  0 

y  0  0  a33  j 


where 


where 


where 


and 


where 


An  =  [Aij]  for 


A 


i  =  1,  •  •  •  ,5,  j  =  1, . . .  ,5, 

0.1,  i=j 
0.0  otherwise 


A22  =  [Ajj]  for 


A 


i  =  1, . . . ,  9,  j  =  1, . . . ,  9, 

1.0,  i  =  j 
0.0  otherwise 


A33  =  [Aij]  for  i  =  1, . . . ,  14,  j  =  1, . . . ,  14, 

.  _  /  2.0,  i=j 

A  ij  —  \ 

0.0  otherwise 


Substituting  the  matrices  T',  A,  and  a  into  Equation  (3.45)  gives  the  failure  time 
distribution  in  the  transform  space.  Application  of  the  one-dimensional  Laplace 
transform  inversion  algorithm  of  Abate  and  Whitt  [1]  is  used  to  arrive  at  the  cumu¬ 
lative  distribution  values  for  selected  points  in  time.  The  results  of  the  failure  time 
distribution  are  compared  to  results  obtained  from  a  simulation  model  in  Table  4.2. 
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Table  4.2  Cumulative  probability  values  for  the  brake  pad  example. 


t 

Simulated 

Analytical 

t 

Simulated 

Analytical 

24.005 

0.009900 

0.009974 

29.855 

0.698370 

0.698667 

24.155 

0.012130 

0.012492 

30.005 

0.720320 

0.720465 

24.305 

0.015210 

0.015505 

30.155 

0.740820 

0.741360 

24.455 

0.019080 

0.019079 

30.305 

0.760130 

0.761315 

24.605 

0.023290 

0.023282 

30.455 

0.779130 

0.780300 

24.755 

0.028480 

0.028184 

30.605 

0.796970 

0.798299 

24.905 

0.034180 

0.033854 

30.755 

0.814180 

0.815301 

25.055 

0.040840 

0.040362 

30.905 

0.830210 

0.831305 

25.205 

0.047980 

0.047774 

31.055 

0.845830 

0.846318 

25.355 

0.055930 

0.056155 

31.205 

0.860470 

0.860351 

25.505 

0.065210 

0.065562 

31.355 

0.873860 

0.873425 

25.655 

0.075500 

0.076048 

31.505 

0.886670 

0.885565 

25.805 

0.086910 

0.087657 

31.655 

0.897460 

0.896800 

25.955 

0.099610 

0.100424 

31.805 

0.908100 

0.907164 

26.105 

0.113240 

0.114374 

31.955 

0.917590 

0.916693 

26.255 

0.127900 

0.129521 

32.105 

0.926270 

0.925428 

26.405 

0.144520 

0.145867 

32.255 

0.934120 

0.933408 

26.555 

0.161680 

0.163401 

32.405 

0.941310 

0.940677 

26.705 

0.181090 

0.182100 

32.555 

0.948080 

0.947278 

27.155 

0.244520 

0.244751 

32.705 

0.953890 

0.953254 

27.305 

0.266340 

0.267615 

32.855 

0.959580 

0.958649 

27.455 

0.290150 

0.291337 

33.005 

0.964200 

0.963504 

27.605 

0.315800 

0.315824 

33.155 

0.968540 

0.967862 

27.755 

0.341430 

0.340974 

33.305 

0.972330 

0.971761 

27.905 

0.366960 

0.366677 

33.455 

0.975700 

0.975241 

28.055 

0.393380 

0.392819 

33.605 

0.978610 

0.978338 

28.205 

0.420010 

0.419282 

33.755 

0.981160 

0.981087 

28.355 

0.445410 

0.445945 

33.905 

0.983770 

0.983520 

28.505 

0.472630 

0.472687 

34.055 

0.985840 

0.985669 

28.655 

0.500000 

0.499389 

34.205 

0.987830 

0.987562 

28.805 

0.526700 

0.525934 

34.355 

0.989590 

0.989224 

28.955 

0.552120 

0.552208 

34.505 

0.990880 

0.990681 

29.105 

0.578320 

0.578105 

34.655 

0.992160 

0.991955 

29.255 

0.603280 

0.603522 

34.805 

0.993230 

0.993065 

29.405 

0.627850 

0.628368 

34.955 

0.994270 

0.994031 

29.555 

0.651890 

0.652556 

35.105 

0.995220 

0.994870 

29.705 

0.675710 

0.676011 

35.555 

0.996890 

0.996760 
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Based  on  this  comparison,  the  technique  developed  in  this  thesis  provides  sim¬ 
ilar  results  to  those  computed  via  simulation.  In  fact,  the  maximum  absolute  devia¬ 
tion  (MAD)  in  probability  for  this  example  is  0.001721.  Figure  4.1  gives  a  graphical 
comparison  of  the  cumulative  distribution  values  at  various  time  points. 


Figure  4.1  Simulated  versus  analytical  cumulative  distribution  functions  for  the  brake 
pad  example. 

Substituting  the  matrices  \Eb  A,  and  a  into  Equation  (3.46)  gives  the  moments 
of  failure  time  in  the  transform  space.  Application  of  the  one-dimensional  Laplace 
transform  inversion  algorithm  of  Abate  and  Whitt  [1]  is  used  to  obtain  the  first 
and  second  moments  of  failure  time,  which  are  compared  to  results  obtained  via 
simulation  in  Table  4.3. 
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Table  4.3  Lower  moments  of  failure  time  for  the  brake  pad  example. 


Measure 

Simulated 

Analytical 

E[TX] 

28.7615 

28.7522 

e  ra 

832.2501 

831.8517 

4-2  Crack  Propagation  in  a  Turbine  Blade 

In  this  example,  the  propagation  of  a  crack  in  a  turbine  blade  of  an  aircraft 
engine  is  considered.  Assume  the  turbine  blade  begins  operation  in  perfect  working 
order  but  develops  a  crack  over  time  due  to  normal  wear. 

Let  X(t)  denote  the  length  of  the  crack  at  time  t.  It  is  assumed  that  the  crack 
grows  at  a  linear  rate  determined  exclusively  by  the  current  state  of  the  random 
environment  to  which  the  turbine  blade  is  exposed.  This  random  environment  can 
be  characterized  as  a  four-state,  semi-Markov  process,  {Z(t)  :  t  >  0},  whose  four 
states  are  defined  by  different  ranges  of  turbine  rotation  speed.  Moreover,  the  rates  at 
which  the  different  states  cause  the  crack  to  propagate  are  known  based  on  previous 
testing.  The  turbine  blade  fails  when  the  crack  reaches  length  x  =  20.0.  The  state 
descriptions,  state  holding  time  distributions,  and  crack  growth  rates  are  listed  in 
Table  4.4.  Additionally,  the  transition  probability  matrix  of  the  embedded  DTMC 
is  given  as 

(  0.0  0.6  0.2  0.2 

0.5  0.0  0.4  0.1 

P  = 

0.2  0.4  0.0  0.4 
v  0.1  0.3  0.6  0.0 

Table  4.4  Description  of  state  space  and  crack  growth  rates  for  the  turbine  blade  ex¬ 
ample. 


State 

Rotation  Speed  (rpm) 

Holding  Time  Distribution 

Crack  Growth  Rate 

1 

1000  through  4999 

Beta(3,5) 

0.17 

2 

5000  through  8999 

Beta(2,4) 

0.43 

3 

9000  through  12999 

Weibull(2,5) 

0.75 

4 

13000  through  15999 

Weibull(3,6) 

1.29 
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From  Table  4.4,  the  matrix  of  degradation  rates  is 


R.£>  — 


(  0.17  0.00  0.00  0.00  N 

0.00  0.43  0.00  0.00 

0.00  0.00  0.75  0.00 

v  0.00  0.00  0.00  1.29  y 


It  is  assumed  that  the  environment  process  begins  in  the  first  state;  therefore,  the 
initial  probability  vector  is 

a=  (  1  0  0  0  )  • 

After  observing  the  system  over  r  =  10000  time  units,  the  T,  matrices  and  T" 
vectors  of  the  phase-type  representations  of  the  four  holding  time  distributions  are 
calculated  as 

ti  =  [<y  f°r  *  =  i,  •  •  • ,  6,j  =  i, . . . ,  6, 


where 


-15.6780, 

15.3140, 

15.6780, 

0 


i  =  j 

i  =  l,j  =  2 

otherwise 


3,. ..,6 


5 


T,  =  [7 i\  for  6, 


where 


7  i  =  < 


0.3638,  i 
0,  i 
15.6780,  i 


1 

2,. ..,5 
6 


T2  =  [5ij]  for  i  =  1, . . . ,  4,  j  =  1, . . . ,  4, 
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where 


where 


where 


where 


where 


-11.6220,  i  =  j 
11.0970,  i  =  1,  j  =  2 
11.6220,  i  —  j  —  1,  j  —  ’ 

0  otherwise 

T2  =  [7*]  for  i  —  1, ...  ,4, 

0.5254,  i  =  1 
7 .  =  <  0,  i  =  2,3, 

11.6220,  i  =  4 

T3  =  [<5ij]  for  i  =  1, . . . ,  4,  j  =  1, . . . ,  4, 

-9.8405,  i  =  j 
9.5442,  i  =  1,  j  =  2 
9.8405,  i  =  j  - 1,  j  =  3,4 
0  otherwise 

T3  =  [7 i]  for  i  —  1, . . . ,  4, 

0.2964,  i  =  1 

7i  =  ■*  0,  i  —  2, ...  ,3  , 

9.8405,  i  =  4 

T4  =  [Sij]  for  i  =  1, . . . ,  8,  j  =  1, . . . ,  8, 

— 16.1640,  i  =  j 
16.0280,  i  =  1,  j  =  2 

16.1640,  /  ./•  1../  3 . 8 

0  otherwise 
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and 


T?4  =  [^j]  for  8, 


where 


7  *  =  < 


0.1357,  i  =  1 
0,  i  —  2, 
16.1640,  i  =  8 


•,7 


Furthermore,  the  observed  transition  rates  among  the  four  environment  states  are 
presented  as  the  elements  of  Q,  where 


-2.6646 

1.5965 

0.5313 

0.5368 

1.5020 

-3.0075 

1.1991 

0.3064 

0.5060 

1.0083 

-2.5170 

1.0027 

0.2033 

0.6087 

1.2234 

-2.0354 

Using  the  above  phase-type  representations  and  the  qVJ  elements  of  Q,  the  26  x  26 
generator  matrix  of  the  newly  formed  CTMC  is 


(  thn  ^12  ^13  ^14  ^ 

^21  ^22  ^23  ^24 

^  = 

^31  Vj/i2  Vi/,.,  Vl/i; 

\  ^41  ^42  ^43  ^44  ) 

where 


^12  =  [Aj\  for  i  =  1,  •  •  • ,  7,  j  =  1, . . . ,  5, 


where 


4>ij 


15965.0,  i  =  7,  j  =  1 

1 

0  otherwise 
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^13  =  [Aj\  for  i  =  1,  •  •  • ,  7,  j  =  1,  •  •  • ,  5, 


where 


5312.6,  i  =  7,j  =  l 
ij) ij  =  <  ,  and 

I  0  otherwise 

^14  =  [Aj\  for  i  =  1,  -  ■ .  ,7,  j  —  1, ...  ,9, 


where 


— 


5368.4,  i  =  7,  j  =  1 
0  otherwise 


^2i  =  [Aj\  for  i  =  1,  •  ■  • ,  5,  j  =  1, . . . ,  7, 


where 


Aj 


15020.0,  i  =  5,  j  =  1 
0  otherwise 


T2 

0  -30075.0 


^23  =  [-0ij]  for  i  =  1,  •  •  • ,  5,  j  =  1, . . . ,  5, 


where 


11991.0,  i  =  5,  j  =  1 

'/>*={ 

|  0  otherwise 

^24  =  [Aj]  for  *  =  1,  •  ■  • ,  5,  j  =  1, . . . ,  9, 


where 


Aj  = 


3063.6,  i  =  5,  j  =  1 
0  otherwise 


^3i  =  [Aj\  for  *  =  1,  ■  •  ■ ,  5,  j  =  1, . . . ,  7, 
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where 


5059.6,  i  =  5,  j  =  1 
0  otherwise 


^32  =  [Aj\  for  i  —  1, ...  ,5,  j  —  1, ...  ,5, 


where 


^ 'ij  ~ 


10083.0,  i  =  5,  j  =  1 
0  otherwise 


T3  T° 

0  -25170.0 


where 


where 


where 


'I' .34  =  [Aj]  for  i  =  1,  ■  •  • ,  5,  j  =  1, . . . ,  9, 


^ ij  ~~ 


10027.0,  i  =  5,  j  =  1 
0  otherwise 


^41  =  [vy  for  i  =  1, . . . ,  9,  j  =  1, . . . ,  7, 


y  = 


2033.2,  i  =  9,  j  =  1 
0  otherwise 


^42  =  [Aj\  for  i  =  1, .  •  ■ ,  9,  j  =  1,  •  •  • ,  5, 


6086.7,  i  =  9,  j  =  1 
0  otherwise 


^43  =  [y]  for  i  =  1, . . . ,  9,  j  =  1, . . . ,  5, 


where 


y  = 


12234.0,  i  =  9,  j  =  1 
0  otherwise 
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and 


^44  = 


t4 

0  -20354.0 


The  26  x  26  expanded  degradation  rate  matrix  is 


where 


where 


where 


where 


and 


A 


(  An  0  0  0  ^ 

0  A22  0  0 

0  0  A33  0 

^  0  0  0  A44  j 


An  =  [Ajj]  for  i  =  1, . . . ,  7,  j  =  1, 


A  ij  — 


1.7,  i=j 
0  otherwise 


A22  =  [Aij]  for  5 ,j  =  1, 


A  ij  — 


4.3,  i  =  j 

1 

0  otherwise 
A33  =  [Kj\  for  i  =  1,  •  •  • ,  5,  j  =  1, 


A«j  — 


7.5,  i  =  j 
0  otherwise 


A44  =  [Ay]  for  i  =  1, . . . ,  9,  j  =  1, 


where 


A, 


12.9,  i=j 
0  otherwise 
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Substituting  the  matrices  \h,  A,  and  a  into  Equation  (3.45)  gives  the  solution  of  the 
failure  time  distribution  in  the  transform  space.  Application  of  the  one-dimensional 
Laplace  transform  inversion  algorithm  of  Abate  and  Whitt  [1]  is  again  used  to  arrive 
at  the  cumulative  distribution  values  for  selected  points  in  time.  The  results  of  the 
failure  time  distribution  are  compared  to  those  obtained  from  a  simulation  model  in 
Table  4.5.  Based  on  this  comparison,  the  technique  developed  in  this  thesis  provides 
similar  results  to  those  computed  via  simulation.  In  fact,  the  MAD  in  probability  for 
this  example  is  0.003297.  Figure  4.2  gives  a  graphical  comparison  of  the  cumulative 
distribution  values  at  various  time  points.  Substituting  the  matrices  fL,  A,  and 


f 


Figure  4.2  Simulated  versus  analytical  cumulative  distribution  functions  for  the  turbine 
blade  example. 

a  into  Equation  (3.46)  gives  the  moments  of  failure  time  in  the  transform  space. 
Application  of  the  one-dimensional  Laplace  transform  inversion  algorithm  of  Abate 


4-17 


Tabic  4.5  Cumulative  probability  values  for  the  turbine  blade  example. 


t 

Simulated 

Analytical 

t 

Simulated 

Analytical 

25.605 

0.016640 

0.016979 

32.005 

0.765210 

0.767488 

25.805 

0.021170 

0.021731 

32.205 

0.788620 

0.791122 

26.005 

0.026420 

0.027490 

32.405 

0.811180 

0.813169 

26.205 

0.032860 

0.034391 

32.605 

0.831400 

0.833613 

26.405 

0.040910 

0.042564 

32.805 

0.850860 

0.852463 

26.605 

0.050410 

0.052139 

33.005 

0.869020 

0.869743 

26.805 

0.061490 

0.063234 

33.205 

0.884060 

0.885496 

27.005 

0.074020 

0.075958 

33.405 

0.898080 

0.899777 

27.205 

0.088680 

0.090403 

33.605 

0.910950 

0.912653 

27.405 

0.104930 

0.106641 

33.805 

0.923320 

0.924198 

27.605 

0.123250 

0.124720 

34.005 

0.933560 

0.934495 

27.805 

0.143160 

0.144660 

34.205 

0.942700 

0.943631 

28.005 

0.164710 

0.166455 

34.405 

0.951060 

0.951694 

28.205 

0.187980 

0.190066 

34.605 

0.958330 

0.958773 

28.405 

0.213170 

0.215423 

34.805 

0.964630 

0.964957 

28.605 

0.240960 

0.242425 

35.005 

0.970190 

0.970332 

28.805 

0.269260 

0.270942 

35.205 

0.974880 

0.974981 

29.005 

0.299300 

0.300818 

35.405 

0.978940 

0.978982 

29.205 

0.330490 

0.331868 

35.605 

0.982560 

0.982410 

29.405 

0.362360 

0.363891 

35.805 

0.985580 

0.985332 

29.605 

0.395080 

0.396667 

36.005 

0.987980 

0.987811 

29.805 

0.428160 

0.429965 

36.205 

0.990040 

0.989906 

30.005 

0.461350 

0.463548 

36.405 

0.992030 

0.991668 

30.205 

0.493880 

0.497177 

36.605 

0.993350 

0.993143 

30.405 

0.527590 

0.530617 

36.805 

0.994600 

0.994374 

30.605 

0.561840 

0.563640 

37.005 

0.995600 

0.995395 

30.805 

0.594630 

0.596032 

37.205 

0.996300 

0.996239 

31.005 

0.626510 

0.627594 

37.405 

0.996860 

0.996934 

31.205 

0.656950 

0.658147 

37.605 

0.997450 

0.997503 

31.405 

0.685580 

0.687534 

37.805 

0.997990 

0.997966 

31.605 

0.713420 

0.715622 

38.005 

0.998490 

0.998342 

31.805 

0.740940 

0.742302 

38.205 

0.998750 

0.998644 

and  Whitt  [1]  is  used  to  obtain  the  first  and  second  moments  of  failure  time,  which 
are  compared  to  results  obtained  via  simulation  in  Table  4.6. 
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Tabic  4.6  Lower  moments  of  failure  time  for  the  turbine  blade  example. 


Measure 

Simulated 

Analytical 

E[T%\ 

30.3253 

30.3046 

SKI 

925.2438 

924.0687 

4-3  Chemical  Decomposition  of  an  Automotive  Coating 

In  this  final  example,  consider  an  automotive  coating  that  is  subject  to  weath¬ 
ering  by  the  outdoor  environment.  Environmental  effects  such  as  temperature  and 
solar  radiation  cause  chemical  decomposition  of  the  coating,  which  can  be  measured 
in  terms  of  gloss  loss  and/or  color  change  [7].  When  the  cumulative  decomposition  of 
the  coating  reaches  a  certain  threshold  x,  the  coating  is  said  to  have  reached  failure. 

Assume  the  chemical  decomposition  (degradation)  rate  is  linear  and  depends 
strictly  on  the  random  state  of  the  outdoor  weather  environment  to  which  the  coating 
is  exposed.  When  the  cumulative  chemical  decomposition  of  the  coating  reaches  the 
value  x  =  5.0,  the  coating  is  said  to  have  failed.  The  weather  environment,  whose 
next  state  depends  only  upon  the  current  state,  can  be  characterized  as  a  semi- 
Markov  process  and  consists  of  five  distinct  states.  These  states,  their  definitions,  the 
probability  distributions  of  their  holding  times,  and  their  associated  decomposition 
rates  are  displayed  in  Table  4.7.  Additionally,  the  transition  probability  matrix  of 


Table  4.7  Description  of  state  space  and  decomposition  rates  for  the  coating  example. 


State 

Sky 

Temp 

Holding  Times 

Decomposition  Rate 

1 

Cloudy 

<  32  deg  F 

Weibull(3,5) 

0.46 

2 

Cloudy 

>  32  deg  F 

Beta(2,5) 

0.82 

3 

Sunny 

<  32  deg  F 

Weibull(4,6) 

1.10 

4 

Sunny 

>  32  deg  F 

Beta(6,3) 

1.34 

5 

Rain 

>  32  deg  F 

Gamma(.5,.l) 

1.98 
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the  embedded  DTMC  is  given  as 


0.5  0.2  0.2  0.1 
0.0  0.3  0.1  0.1 
0.4  0.0  0.2  0.2 
0.1  0.2  0.0  0.6 
0.1  0.1  0.6  0.0 

From  Table  4.7,  the  matrix  of  degradation  rates  is 

(  0.46  0.00  0.00  0.00 
0.00  0.82  0.00  0.00 
R D  =  0.00  0.00  1.10  0.00 

0.00  0.00  0.00  1.34 

v  0.00  0.00  0.00  0.00 

Again,  it  is  assumed  that  the  environment  process  begins  in  the  first  state;  therefore, 
the  initial  probability  vector  is 

a  =  (  1  0  0  0  o). 

After  observing  the  system  over  r  =  40000  time  units,  the  T,  matrices  and  T ° 
vectors  of  the  phase-type  representations  of  the  five  holding  time  distributions  are 
calculated  as 

Ti  =  [Sij]  for  i  =  1, . . . ,  6,  j  =  1, . . . ,  6, 

where 

-15.694,  i  =  j 
15.328,  i  =  1,  j  =  2 
15.694,  i=j- l,j  =  3,...,6  ’ 

0  otherwise 
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T,  =  [7 i]  for  6, 


where 


where 


where 


0.36645,  i 


7 i  =  < 


0,  i 


15.694,  i 


1 

2,. ..,5 
6 


T2  =  [<y  for  i  =  1, . . . ,  4,  j  =  1, . . . ,  4, 


-9.8485, 

9.602, 

9.8485, 

0 


i  =  j 

i  =  l,j  =  2 
i=j  -1,3 

otherwise 


3,4 


5 


T"  =  [7 i]  for  4, 


0.24652,  i  =  1 
7 *  =  <  0,  i  =  2, . . . ,  3 

9.8485,  i  =  4 

T3  =  [5*j]  for  i  =  l,...,8,j  =  1,. 


where 


where 


f  -19.704,  i  =  j 


19.444, 

19.704, 

0 


*  =  1,  j  =  2 

*  =  J  -  !,  J  =  3,  •  •  • ,  8 
otherwise 


T 


O 

3 


7 i]  for  *  =  8, 


5 


0.25978,  i 


7*  =  < 


0,  i 


19.704,  i 


1 


2,. ..,7 


5 


8 
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T4  =  [Sij]  for  i  =  1, . . . ,  8,  j  =  1, . . . ,  8, 


where 


where 


where 


and 


where 


-16.122,  i—j 


16.01, 

16.122, 

0 


*  =  1,  j  =  2 

*  =  j  —  1,  j  =  3,  •  •  • ,  8 

otherwise 


T?4  =  [^j]  for  8, 


5 


%  =  < 


0.11148,  i  =  1 
0,  i  =  2,...,  7 
16.122,  i  =  8 


T5  =  for  i  =  1, . . . ,  2,  j  =  1, . . . ,  2, 


-22.404,  i  =  j 


= 


0.43763, 

-3.6311, 

0 


*  —  1,4—2 
i  =  2,j  =  2 
otherwise 


5 


Tg  =  N  for  *  =  1,2, 

(  21.967,  i  =  1 
h  =  < 

3.6311,  i  =  2 
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Furthermore,  the  observed  transition  rates  among  the  five  environment  states  are 
presented  as  the  elements  of  Q,  where 


-2.6675 

1.3358 

0.5379 

0.5336 

0.2603 

1.2608 

-2.5092 

0.7500 

0.2502 

0.2483 

0.5079 

0.9927 

-2.4917 

0.5004 

0.4907 

0.2001 

0.2033 

0.4167 

-2.0275 

1.2074 

4.0752 

1.9409 

1.9833 

11.9950 

-19.9940 

Using  the  above  phase- type  representations,  the  ql3  elements  of  Q,  and  M  =  10000, 
the  33  x  33  generator  matrix  of  the  newly  formed  CTMC  is 


where 


/ 

^11 

*12 

* 

CO 

*14 

*15  ^ 

^21 

*22 

*23 

*24 

*25 

^31 

*32 

•6 

CO 

CO 

*34 

*35 

^41 

*42 

*43 

*44 

*45 

V 

*51 

*52 

*53 

*54 

*55  J 

= 


-26675.0  j 


5 


^12  =  [Aj\  for  *  =  1,  ■  ■  ■ ,  7,  j  =  1, . . . ,  5, 


where 


^ ij 


13358.0,  i  =  7,  j  =  1 

1 

0  otherwise 


^13  =  [Aj\  for  i  =  1,  •  •  • ,  7,  j  =  1, . . . ,  9, 


where 


'Ipij 


5378.7,  i  =  7,  j  =  1 
0  otherwise 


,  and 
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where 


where 


^14  =  [Aj\  for  i  =  1, . . . ,  7,  j  =  1, . . . ,  9, 


Aj  = 


5335.6,  i  =  7 ,j  =  1 
0  otherwise 


^15  =  [Aj]  for  i  =  1,  -  ■  ■ ,  7,  j  =  1, . . . ,  3, 


Aj  — 


2602.6,  i  —  7,  j  —  1 
0  otherwise 


^21  =  [Aj]  for  i  =  1,  •  ■  • ,  5,  j  =  1, . . . ,  7, 


where 


Aj 


12608.0,  i  =  5,  j  =  1 
0  otherwise 


where 


where 


0  -25092.0 


^23  =  [Aj]  for  i  =  1,  •  •  • ,  5,  j  =  1, . 


Aj  ~ 


7499.9,  i  =  5,  j  =  1 
0  otherwise 


^24  =  [V’ij]  for  i  =  1, .  •  ■ ,  5,  j  =  1, . . . ,  9, 


Aj  = 


2502.0,  i  =  5,  j  =  1 
0  otherwise 


^25  =  [V’ij]  for  i  =  1, . . . ,  5,  j  =  1, . . . ,  3, 


where 


2482.5,  i  =  5,j  =  l 

Aj  —  \  ) 

0  otherwise 


4-24 


where 


^3i  =  [Aj]  for  i  —  1, . . ,  ,9,  j  —  1, ...  ,7, 


Aj  = 


5079.3,  i  =  9,  j  =  1 
0  otherwise 


^32  =  [Aj]  for  i  =  1,  •  •  • ,  9,  j  =  1,  •  •  • ,  5, 


where 


Aj  — 


9926.7,  i  =  9,  j  =  1 
0  otherwise 


T3  T° 

0  -24917.0 


where 


where 


^34  =  [Aj]  for  i  =  1,  •  •  • ,  9,  j  =  1, . . . ,  9, 


V’ij  — 


5004.1,  i  =  9,  j  =  1 
0  otherwise 


^35  =  [Aj]  for  i  =  1,  •  •  • ,  9,  j  =  1,  •  •  • ,  3, 


= 


4907.3,  i  =  9,  j  =  1 
0  otherwise 


^4i  =  [t/'ij]  for  i  =  1, . . . ,  9,  j  =  1,  • .  ■ ,  7, 


where 


2000.5,  i  =  9,  j  =  1 
Aj  =  \  j 

0  otherwise 


^42  =  [V’ij]  for  i  =  1,  •  •  • ,  9,  j  =  1, . . . ,  5, 
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where 


Aj  = 


2033.4,  %  =  9,  j  =  1 
0  otherwise 


^43  =  [Aj]  for  i  —  1, ...  ,9,  j  —  1, ...  ,9, 


where 


Aj  — 


4167.4,  i  =  9,  j  =  1 
0  otherwise 


where 


where 


where 


T4  T° 

0  -20275.0 


vh 45  =  [Aj\  for  i  =  1,  •  •  • ,  9, 3  =  1,  • 


Aj 


12074.0,  i  =  9,  j  =  1 
0  otherwise 


^5i  =  [Aj]  for  i  =  1,  •  •  • ,  3,  j  =  1, . . . ,  7, 


Aj  — 


40752.0,  i  =  3,  j  =  1 
0  otherwise 


^52  =  [''/y  for  i  =  1, . . . ,  3,  j  =  1, . . . ,  5, 


Aj 


19409.0,  i  =  3,  j  =  1 
0  otherwise 


^5.3  =  [y  ]  for  i  =  1, . . . ,  3,  j  =  1, . . . ,  9, 


where 


19833.0,  i  =  3,  j  =  1 
ipij  =  <  ,  and 

0  otherwise 


^54  =  [Aj\  for  i  =  1,  - .  ■ ,  3,  j  =  1, . . . ,  9, 
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where 


Aj  = 


119950.0000,  i  =  3,  j  =  1 
0  otherwise 


0  -199940.0 


The  33  x  33  expanded  degradation  rate  matrix  is 


An  0  0  0  0 

0  A22  0  0  0 

0  0  A33  0  0 

0  0  0  A44  0 

0  0  0  0  A5; 


where 


where 


where 


An  =  [Kj\  for  *  =  1,  •  •  • ,  7,  j  =  ljn  ■ ,  7, 


Ay  — 


2.3,  i  =  j 
0  otherwise 


A22  =  [Aij]  for  i  =  1, . . . ,  5,  j  =  1, . . . ,  5, 


A*j  — 


4.1,  i=j 
0  otherwise 


A 33  =  [Xij]  for  i  =  1, . r9,  j  =  1, . . . ,  9, 


where 


A  i? 


5.5,  i  =  j 
0  otherwise 
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and 


A44  =  [Xij]  for  i  =  1, . . . ,  9,  j  =  1, . . . ,  9, 


where 


and 


where 


A 


6.7,  i=j 
0  otherwise 


A55  =  [Aij]  for  i  =  1, . . . ,  3,  j  =  1, . . . ,  3, 


9.9,  %  =  j 
0  otherwise 


Applying  th,  A,  and  a  to  Equation  (3.45)  gives  the  solution  of  the  failure  time  distri¬ 
bution  in  the  transform  space.  The  analytical  results  of  the  failure  time  distribution 
in  the  time  domain  are  compared  to  results  derived  via  simulation  in  Table  4.8. 
Based  on  this  comparison,  the  technique  developed  in  this  thesis  provides  similar 
results  to  those  computed  via  simulation.  In  fact,  the  MAD  in  probability  for  this 
example  is  0.004181.  Figure  4.3  gives  a  graphical  comparison  of  the  cumulative  dis¬ 
tribution  values  at  various  time  points. 
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Table  4.8  Cumulative  probability  values  for  the  coating  example. 


t 

Simulated 

Analytical 

t 

Simulated 

Analytical 

4.155 

0.021700 

0.022164 

5.640 

0.726490 

0.723197 

4.210 

0.029530 

0.029350 

5.695 

0.754040 

0.750956 

4.265 

0.038490 

0.038055 

5.750 

0.780140 

0.777106 

4.320 

0.049100 

0.048431 

5.805 

0.804230 

0.801577 

4.375 

0.061180 

0.060610 

5.860 

0.827300 

0.824326 

4.430 

0.075030 

0.074701 

5.915 

0.847530 

0.845333 

4.485 

0.091160 

0.090782 

5.970 

0.866250 

0.864603 

4.540 

0.109060 

0.108900 

6.025 

0.883390 

0.882159 

4.595 

0.128770 

0.129070 

6.080 

0.898950 

0.898048 

4.650 

0.150780 

0.151274 

6.135 

0.912910 

0.912330 

4.705 

0.174850 

0.175464 

6.190 

0.925670 

0.925079 

4.760 

0.201460 

0.201560 

6.245 

0.937570 

0.936383 

4.815 

0.228750 

0.229453 

6.300 

0.947580 

0.946335 

4.870 

0.259000 

0.259001 

6.355 

0.955760 

0.955035 

4.925 

0.291560 

0.290033 

6.410 

0.963040 

0.962588 

4.980 

0.323730 

0.322352 

6.465 

0.968900 

0.969095 

5.035 

0.356320 

0.355736 

6.520 

0.974440 

0.974662 

5.090 

0.390980 

0.389948 

6.575 

0.978930 

0.979388 

5.145 

0.425950 

0.424735 

6.630 

0.983290 

0.983370 

5.200 

0.461860 

0.459840 

6.685 

0.986990 

0.986698 

5.255 

0.496930 

0.495004 

6.740 

0.989770 

0.989457 

5.310 

0.532100 

0.529974 

6.795 

0.991880 

0.991726 

5.365 

0.567740 

0.564506 

6.850 

0.993720 

0.993576 

5.420 

0.601350 

0.598370 

6.905 

0.995140 

0.995071 

5.475 

0.634620 

0.631353 

6.960 

0.996270 

0.996268 

5.530 

0.667410 

0.663261 

7.015 

0.997180 

0.997218 

5.585 

0.697890 

0.693925 

7.070 

0.997820 

0.997964 
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Figure  4.3  Simulated  versus  analytical  cumulative  distribution  functions  for  the  coating 
example. 

Substituting  the  matrices  \Eb  A,  and  a  into  Equation  (3.46)  gives  the  moments 
of  failure  time  in  the  transform  space.  Application  of  the  one-dimensional  Laplace 
transform  inversion  algorithm  of  Abate  and  Whitt  [1]  is  used  to  obtain  the  first 
and  second  moments  of  failure  time,  which  are  compared  to  results  obtained  via 
simulation  in  Table  4.9. 

Table  4.9  Lower  moments  of  failure  time  for  the  coating  example. 


Measure 

Simulated 

Analytical 

E[T,\ 

5.2830 

5.2867 

m a 

28.2738 

28.3145 

The  examples  in  this  chapter  show  that  the  failure  time  distribution  of  a  system 
experiencing  wear  due  to  an  environment  characterized  as  a  semi-Markov  process  can 
be  approximated  by  simply  observing  the  environment  over  time.  In  all  three  nu¬ 
merical  examples,  the  maximum  absolute  deviation  (MAD)  in  probability  was  less 
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than  or  equal  to  0.0042.  Considering  the  inherent  inaccuracies  due  to  the  use  of 
phase-type  approximations,  the  results  are  favorable.  Assuming  that  the  wear  rates 
imposed  on  the  system  by  each  state  are  known,  lifetime  estimation  can  be  done 
without  monitoring  the  cumulative  degradation  level  during  the  system’s  operating 
lifetime.  This  environment-based  method  is  particularly  advantageous  over  tradi¬ 
tional  degradation-based  reliability  techniques  when  the  operating  environment  of 
the  system  is  known  or  observable  but  degradation  measurements  are  difficult  or 
impossible  to  obtain. 
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5.  Conclusions  and  Future  Research 

This  thesis  presented  a  reliability  analysis  technique  for  a  system  that  expe¬ 
riences  wear  under  the  influence  of  a  dynamic  environment  process.  In  particular, 
this  effort  provided  a  means  to  compute  the  failure  time  distribution  and  moments 
of  failure  time  for  a  single-unit  system  whose  wear  rate  depends  on  the  state  of 
its  stochastic  environment  when  the  environment  is  characterized  as  a  semi-Markov 
process.  The  numerical  examples  showed  that  a  lifetime  prediction  can  be  made  for 
a  degrading  system  simply  by  observing  the  system’s  operating  environment  over 
time. 

The  first  step  in  this  research  effort  was  to  review  the  two  main  approaches  to 
system  lifetime  estimation:  classical  failure  time  analysis  and  degradation-based  reli¬ 
ability  analysis.  After  a  more  through  study  of  the  current  literature  on  degradation- 
based  reliability  techniques,  it  was  observed  that  a  common  drawback  of  these  ap¬ 
proaches  is  the  lack  of  consideration  of  the  effects  of  the  operating  environment  on  the 
performance  of  the  system.  Since  many  systems  operate  under  complex  and  chang¬ 
ing  conditions,  a  study  of  the  effects  of  the  environment  on  a  system  is  necessary 
to  obtain  a  realistic  estimation  of  the  system’s  reliability.  One  environment-based 
model  produced  good  results,  as  compared  to  those  from  a  simulation  model,  when 
the  state-dependent  wear  rates  were  known  and  the  environment  was  characterized 
as  a  finite  state  continuous-time  Markov  chain  (CTMC)  [16].  However,  this  model 
assumed  that  the  time  spent  in  each  of  the  environment  states  is  an  exponentially 
distributed  random  variable.  Under  realistic  conditions,  the  state  holding  time  dis¬ 
tributions  may  not  be  exponentially  distributed  or  even  known. 

If  the  evolution  of  the  environment  is  assumed  to  depend  only  upon  the  current 
state,  then  any  such  environment  process  has  an  embedded  discrete-time  Markov 
chain  (DTMC)  and  can  be  characterized  as  a  semi-Markov  process.  In  such  cases, 
knowledge  of  the  holding  times  in  each  state  may  be  gained  solely  by  observation. 
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This  thesis  presented  a  technique  in  which  an  environment  characterized  as  a  semi- 
Markov  process  can  be  re-characterized  as  a  continuous-time  Markov  chain.  This  was 
accomplished  by  approximating  the  state  holding  time  distributions  with  phase-type 
distributions.  An  overview  of  phase-type  approximation  techniques  was  presented  in 
chapter  3.  The  transition  rates  among  the  various  states  of  the  environment  and  the 
phase-type  representations  were  then  used  to  construct  the  infinitesimal  generator 
matrix  of  the  newly  created  CTMC.  This  generator  matrix,  the  appropriately  resized 
matrix  of  wear  rates,  and  the  initial  probability  vector  of  the  environment  process 
were  used  to  compute  the  failure  time  distribution  and  moments  of  failure  time  of 
the  system  in  the  transform  using  the  one-dimensional  technique  found  in  [17]. 

After  the  entire  procedure  was  outlined  in  chapter  3,  numerical  examples  were 
illustrated  in  chapter  4  to  show  the  potential  accuracy  of  failure  time  predictions 
using  the  phase-type  approximation  technique.  In  each  example,  the  environment 
was  observed  for  a  sufficiently  long  time  period  during  which  the  number  of  transi¬ 
tions  among  the  various  states  and  the  state  holding  times  were  observed.  From  the 
transition  observations,  transition  rates  were  estimated  for  each  state  of  the  environ¬ 
ment.  Phase-type  representations  of  the  state  holding  time  distributions,  computed 
using  an  algorithm  implemented  in  MATLAB®,  were  then  used  with  the  transition 
rates  to  construct  the  infinitesimal  generator  matrix  of  the  newly  formed  CTMC. 
The  degradation  rate  matrix  was  expanded  to  account  for  the  number  of  phases  used 
to  approximate  each  state.  The  expanded  degradation  rate  matrix,  the  generator 
matrix,  and  the  initial  probability  vector  of  the  environment  process  were  then  used 
as  inputs  to  Equations  (3.45)  and  (3.46)  to  obtain  the  failure  time  distribution  val¬ 
ues  and  moments  of  failure  time,  respectively,  in  the  transform  space.  A  numerical 
inversion  technique  [1],  implemented  in  MATLAB®,  was  then  used  to  obtain  dis¬ 
tribution  values  and  moments  in  the  time  domain.  These  values  were  compared  to 
those  obtained  via  simulation  to  show  the  relative  accuracy  of  lifetime  estimation.  In 
all  three  examples,  the  maximum  absolute  deviation  (MAD)  in  probability  was  less 
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than  or  equal  to  0.0042.  Considering  that  the  infinitesimal  generator  matrix  used 
in  Equation  (3.45)  was  based  upon  phase- type  approximations  of  the  state  holding 
time  distributions,  these  comparisons  provide  very  favorable  results. 

The  main  contribution  of  this  thesis  effort  is  the  resulting  procedure  for  mak¬ 
ing  lifetime  predictions  of  degrading  systems  based  solely  upon  observations  of  their 
operating  environments.  In  particular,  this  work  generalizes  the  results  in  [16]  and 
[17]  to  the  case  where  state  holding  time  distributions  are  nonexponential  or  even 
unknown.  By  using  phase-type  distributions  to  approximate  state  holding  time 
distributions,  the  method  developed  in  this  thesis  eliminates  the  requirement  for 
exponentially  distributed  holding  times,  allowing  any  environment  that  can  be  char¬ 
acterized  as  a  semi-Markov  process  to  be  converted  into  a  CTMC.  Furthermore,  this 
thesis  expands  the  knowledge  base  of  degradation-based  reliability  analysis  by  pro¬ 
viding  a  technique  for  failure  time  estimation  that  does  not  rely  upon  degradation 
measurements  during  the  lifetime  of  the  system.  This  technique  can  prove  partic¬ 
ularly  useful  in  situations  where  the  system’s  operating  environment  is  observable 
or  known,  but  where  observations  of  the  cumulative  wear  to  the  system  over  time 
are  not  easy  (or  even  possible)  to  obtain.  For  example,  designers  of  components  on 
satellites  or  systems  that  operate  under  deep  water  may  benefit  from  this  technique. 
Such  systems  may  operate  in  predictable  or  observable  environments,  but  restric¬ 
tions  due  to  weight,  cost,  or  environmental  conditions  may  prevent  the  collection  of 
degradation  data. 

Suggested  areas  for  future  research  include  expansion  of  the  phase-type  ap¬ 
proximation  technique  to  those  other  than  the  Coxian  and  Erlang  distributions. 
These  distributions  were  used  due  to  their  simplicity  and  the  fact  that  all  arbitrary 
distributions  can  be  approximated  by  a  Coxian  or  Erlang  distribution  based  on  the 
coefficient  of  variation  of  the  original  distribution.  The  use  of  other  phase-type  distri¬ 
butions  may  increase  the  accuracy  of  the  approximation  when  compared  to  empirical 
or  known  parametric  distributions.  This  improved  approximation  technique  should 
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enhance  the  accuracy  of  the  lifetime  estimation  when  compared  to  simulated  results. 
Furthermore,  modifications  of  this  procedure  that  reduce  the  number  of  phases  used 
to  approximate  each  sojourn  time  distribution  would  greatly  reduce  computational 
effort  due  to  a  dimensionality  reduction  of  the  infinitesimal  generator  matrix  and  as¬ 
sociated  degradation  rate  matrix.  In  their  recent  work  on  phase-type  distributions, 
Osogami  and  Harchol-Balter  [22]  present  a  moment-matching  algorithm  that  nearly 
minimizes  the  number  of  phases  used  to  approximate  an  arbitrary  distribution.  The 
time  available  for  observation  of  the  system’s  operating  environment,  and  an  accept¬ 
able  balance  between  the  accuracy  of  the  phase-type  approximations  and  the  number 
of  phases  used  for  each  state,  will  determine  the  efficacy  of  the  procedures  developed 
in  this  thesis. 
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Appendix  A.  Code  for  Brake  Pad  Example 


%  Plot_SMP_weib .m 

o, 

%  This  MATLAB  program  runs  the  semi-Markov  process  (SMP)  simulation 
(semi_Markov_weib . m)  and  then  the  'simulation  plus  phase-type  approximations 
plus  analytical  solution'  program  SMPweibstates .m.  The  CDF  values  are  then 
compared  in  a  plot  (Figure  2) . 

o, 

o, 

%  Author:  Captain  Christopher  J.  Solo,  M.S.  candidate,  Operations  Research 
%  Air  Force  Institute  of  Technology 

%  Date:  January  29,  2004 

o, 

o, 

clear  all; 

o, 

t=0 . 005 : 0 . 05 : 50 . 0 ;  %  vector  of  time  points  at  which  to  evaluate  both  CDFs 

o, 

semi_Markov_weib;  %  calls  program  to  simulate  semi-Markov  process  (SMP)  and 
obtain 

%  failure  times 

O, 

SMPweibstates;  %  calls  program  to  simulate  SMP  (long  term),  compute  phase-type 
%  approximations,  and  produce  analytical  results  of  failure 

time  distribution 

o. 

%  F  =  cdf  values  from  simulation 

%  FailureTimeProb  =  cdf  values  from  phase-type  analytical  solution 

o, 

%  Produces  plot  of  simulated  failure  time  distribution  vs.  analytical 
%  CDF  derived  via  conversion  of  environment  process  from  SMP  to  CTMC 

o, 

figure  (2) ; 
plot (t , F, ' r ' ) ; 
hold  on; 

plot (t ,  FailureTimeProb,  ' b ' ) ; 
grid  off; 

title (' Failure-time  CDF  values:  Simulated  vs.  Analytical'); 
xlabel ( ' t '  )  ; 
ylabel ( 'P (T  <  t)  '  )  ; 

legend (' Simulated  CDF ',' Analytical  CDF',0); 


A-l 


PROGRAM  semi_Markov_weib . m 

The  purpose  of  this  MATLAB  program  is  to  simulate  a  finite-state  semi-Markov 
process.  The  process  is  simulated  in  order  to  validate  analytical  solutions 
for  the  distribution  and  moments  of  the  lifetime  of  components  whose 
degradation  threshold  is  WearThreshold .  The  program  uses  function  "rando"  in 
order  to  select  the  next  state  after  a  state  transition. 


Orig  Author: 

Date : 
Revised  by: 

Date  : 


Jeffrey  P.  Kharoufeh,  Ph.D.  candidate,  IE  &  OR, 
Penn  State  University 
January  15,  2001 

Captain  Chris  Solo,  M.S.  candidate,  OR, 

Air  Force  Institute  of  Technology 
January  29,  2004 


VARIABLE  DEFINTIONS 


if-k-k-k-k-k'k'k-k-k'k'k-k-k-k'k'k'k'k-k'k-k-k-k-k-k'k'k-k'k'k'k-k-k'k'k'k-k'k'k-k-k-k'k'k'k-k-k'k-k-k-k'k'k'k'k-k'k-k-k-k-k'k'k'k 


%  WearThreshold  =  Degradation  threshold  (signifies  failure) 
%  k  =  an  index  variable 


%  VECTOR/FUNCTION  DEFINITIONS 

o, 

%  Lifetime  =  Vector  of  lifetimes 

%  B  =  Vector  for  the  initial  probability  distribution  of  the  Markov  process 
%  P  =  Probability  transition  matrix 

%  V  =  Vector  of  degradation  rates  (i.e.,  V(i)=  degradation  rate  when  environment 
%  is  in  state  i) . 

%  Z  =  Vector  of  environment  states 


%  The  program  assumes  a  state  space  of  the  form  S={1,2, . . .K} 


semi_Markov_input_weib;  %  Obtain  intialization  parameters 

Z  =  [];  %  Initialize  Z. 

for  k  =  1 : length (Lifetime) 

Z  =  []  ; 

Z (1) =rando (B) ;  %  Initial  state  of  the  environment  at  time  0 

Lifetime(k)  =  weibrnd (K ( Z ( 1 ) , 1 ) , K ( Z ( 1 ) , 2 ) ) ;  %  Time  spent  in  initial  state 
D  =  0.0;  %  At  time  0,  degradation  is  0 

D  =  D  +  Lifetime (k) *V(Z (1) ) ;  %  Add  degradation  to  cumulative  degradation 

i— l  ; 


while  D  <  WearThreshold  %  Do  while  cumulative  degradation  <  WearThreshold 
Z(i+1)  =  rando (P (Z (i) ,:)) ;  %  Use  the  matrix  P  to  determine  next  state 

new_time  =  weibrnd (K ( Z (i+1) , 1 ), K (Z (i+1 ), 2 )) ;  %  Time  spent  in  state  i+1 
D  =  D  +  new_time*V (Z (i+1) ) ;  %  Update  cumulative  degradation 
Lifetime (k)  =  Lifetime (k)  +  new_time;  %  Update  total  lifetime 

i— i+1 ; 

end 

Lifetime (k) =Lifetime (k) - (1/V (Z (i) ))* (D-WearThreshold) ;  %  Subtract  time  after 

%  reaching  threshold 

end 

get_cdf_weib;  %  Computes  empirical  distribution 
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%  PROGRAM  semi_Markov_input_weib .m 

O, 

%  The  purpose  of  this  MATLAB  program  is  to  provide  input  data  to  program 
%  semi_Markov_weib . m  for  simulation  of  a  finite-state  semi -Markov  process 
%  (K  states) .  The  process  is  simulated  in  order  to  validate  analytical 
%  solutions  for  the  distribution  and  moments  of  the  failure  time  of  systems  with 
%  degradation  threshold  L  subject  to  a  SMP  environment.  The  program  uses 
%  function  "rando"  in  order  to  select  the  next  state  after  a  state  transition. 


%  Orig  Author: 

O, 

%  Date: 

%  Revised  by: 

o, 

%  Last  Revised: 


Jeffrey  P.  Kharoufeh,  Ph.D.  candidate,  IE  &  OR, 
Penn  State  University 
January  15,  2001 

Captain  Christopher  Solo,  M.S.,  OR, 

Air  Force  Institute  of  Technology 
January  17,  2004 


%  VARIABLE  DEFINTIONS 


%  WearThreshold  =  degradation  threshold  (system  fails  after  accumulating  L 
%  amount  of  damage) 

o, 

%  VECTOR/FUNCTION  DEFINITIONS 

o, 

%  T  =  Vector  of  failure  time  observations 

%  B  =  Vector  for  the  initial  probability  distribution  of  the  semi-Markov  process 
%  P  =  Probability  transition  matrix 

%  V  =  Vector  of  degradation  rates  (i.e.,  V(i)=  degrad,  rate  when  environment  is 
%  in  state  i) . 


%  Initialize  variable  and  vector  values 


WearThreshold  =3.0;  %  Component  reaches  'failure'  with  this  level  of 

degradation 

N  =  3;  %  number  of  states 


Lifetime  =  zeros  ( 1 ,  100000 ) ;  %  number  of  simulations  (failure  time 

observations ) 

B  =  [1  zeros (1, N-l )] ; 


S  =  [1,2,3];  %  states  of  the  environment 

V  =  [1,10,20];  %  vector  of  degradation  rates 


P  =  [0.0  0.7  0.3;  %  transition  probability  matrix; 

0.6  0.0  0.4; 

0.2  0.8  0.0] ; 


K  =  [4  2;  %  Parameter  matrix — each  row  gives  parameters  of  distribution 

5  3; 

6  4]; 
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%  PROGRAM  get_cdf_weib . m 

%  The  purpose  of  this  MATLAB  program  is  to  construct  an  empirical  distribution 
%  function  based  upon  a  vector  of  observations  of  some  continuous  random 
%  variable.  The  values  of  the  CDF  are  generally  used  for  the  purpose  of 
%  performing  hypothesis  tests  on  the  equality  of  two  nonparametric  distributions. 
%  Other  approaches  for  testing  may  be  the  Cramer  von  Mises  test,  which  is 
%  similar  to  the  K-S  two-sample  test,  but  slightly  more  robust. 


Orig  Author:  Jeffrey  P.  Kharoufeh,  Ph.D.  candidate,  IE  &  OR, 
Penn  State  University 
Date:  June  2000 

Revised  by:  Captain  Chris  Solo,  M.S.  candidate,  OR, 

Air  Force  Institute  of  Technology 
Last  Revised:  January  18,  2004 


%  VECTOR  DEFINITIONS 


%  Lifetime  =  The  finite-dimensional  vector  of  real  observations  (read  from 
%  semi_Markov_weib .m) 

%  t  =  The  vector  of  points  at  which  to  estimate  the  CDF. 

%  F  =  The  vector  of  CDF  estimates.  That  is,  F(i)  =  F(t(i)),  i=l,2,...,n  of 
%  U,2, ...,C}. 


F  =  zeros ( 1 , length (t )) ;  %  This  ensures  that  t  and  F  are  vectors  of  equal  length 

%  Compute  the  cdf  value  at  the  point  tO 

for  i  =  1: length (t) 

F (i)  =  0.0; 

for  j  =  1 : length (Lifetime) 
if  Lifetime ( j ) <=  t (i) 

F(i)  =  F(i)  +  1/length (Lifetime) ; 

else 

F  (i)  =  F  ( i )  ; 

end 

end 

end 

%  WRITE  OUTPUTS  TO  A  TAB  DELIMITED  TEXT  FILE 

fid  =  fopen ( ' F : \AFIT\Thesis\Solo  Thesis\SMP 
s imulation\ re suit s\brake simulated . out '  ,  ' w ' ) ; 
for  b=l : length (t) 

fprintf (fid,  ' %4 . 6f\t%4 . 6f\n' , t (b) , F (b) ) ; 

end 

f close (fid) 
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%  SMPweibstates  .m 

o, 

%  The  purpose  of  this  MATLAB  program  is  to  simulate  a  semi-Markov  process  (SMP) 

%  on  a  state  space  S  and  convert  it  into  a  continuous-time  Markov  chain,  using  a 
%  phase-type  distribution  to  approximate  the  holding  time  distribution  of  each 
%  environment  state.  Each  of  the  environment  states  is  known  to  have  a  Weibull 
%  holding  time  distribution  (each  with  different  parameters.)  The  output 
%  is  the  analytical  solution  for  the  failure  time  distribution  of  a  system 
%  subject  to  the  SMP  environment,  based  on  the  results  of  Kharoufeh  (2003)  and 
%  Kharoufeh  and  Sipe  (2004)  . 

o, 

%  Original  code  for  process. m  and  rando.m  obtained  from  Craig  L.  Zirbel  at 
%  http://www-math.bgsu.edu/~zirbel/ap/  It  appeared  that  the  code  was  based 
%  on  a  discrete  time  Markov  chain,  and  Captain  Steve  Cox  modified  it  to  work 
%  for  a  continuous-time  Markov  chain.  Captain  Chris  Solo  modified  it  to  be 
%  a  SMP  (January  2004) . 

o, 

o, 

%  Author:  Captain  Christopher  J.  Solo,  M.S.  candidate, 

%  Operations  Research,  Air  Force  Institute  of  Technology 

%  Date:  September  13,  2003 

%  Last  Revised:  January  18,  2004 

o, 

%  USER  INPUTS 

%  Tmax  =  time  of  SMP  run 

%  numruns  =  number  of  simulations  (runs) 

%  mu  =  initial  probability  vector  of  environment  states 
%  S  =  state  space  vector 

%  R  =  degradation  rate  vector  of  environment  states 
%  P  =  transition  probability  matrix 

%  K  =  matrix  of  Weibull  parameters  for  environment  states 
%  M  =  large  value  used  in  Big  M  (hot  potato)  method 

o, 

%  Additionally,  rando.m,  phasetypeapprox . m,  PHvalues.m,  invt_lap_Solo . m,  and 
%  failLST.m  are  needed. 

o, 

%  For  many  short  runs,  set  numruns  high  and  Tmax  low.  For  one  long  run, 

%  set  numruns  equal  to  1  and  make  Tmax  high. 

format  long  g; 

%  User  Inputs 

o, 

numruns  =  1;  %  the  number  of  times  (runs)  to  simulate  the  process 

Tmax  =  10000;  %  This  is  the  length  of  the  run  of  the  environment  process 

mu  =  [1,0, 0,0];  %  Semi-Markov  process  (SMP)  starts  in  state  1 

S  =  [1,2,3];  %  There  are  three  states  defined 

R  =  [.1,1,2];  %  Degradation  rate  of  each  state  of  environment 
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P  =  [0.0  0.7  0.3;  %  transition  probability  matrix; 

0.6  0.0  0.4; 

0.2  0.8  0.0] ; 

K  =  [4  2;  %  Parameter  matrix — each  row  gives  parameters  of  different 

distribution 

5  3; 

6  4]; 

^•k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

%  Simulation  of  the  semi -Markov  process  (SMP) 

statetimeTij  =  zeros (length (S) , length  (S) ) ; 
statetimeTi  =  []; 
scale  =  1; 

TimeToFail ( 1 )  =0;  %  start  times  at  0 

%Tall (runnum, 1)  =  0; 

x(l)  =  rando (mu) ;  %  determine  starting  state  (time  0,  not  time  1) 

%xall (runnum, 1)  =  x(l); 
rates(l)  =  R(x(l)); 

%ratesall (runnum, 1)  =  rates  (1); 

i  =  1; 

o, 

while  TimeToFail (i)  <  Tmax, 

x(i+l)  =  rando (P (x (i) ,:)) ;  %  Row  vector  from  P-matrix 

determines  next  transitions 

%xall (runnum, i+1 )  =  x(i+l); 

time(i)  =  weibrnd (K (x ( i ) , 1 ) , K (x ( i) , 2 ) ) ;  %  Generate  Weibull  rv  based  on 

%  current  state 

statetimeTi j (x (i) , x (i+1) )  =  statetimeTi j (x (i) , x (i+1 ) )  +  time(i); 

TimeToFail (i+1)  =  TimeToFail (i)  +  time(i);  %  Add  to  current  time 
%Tall (runnum, i+1 )  =  T(i+1); 

rates  (i  +  1)  =  R(x(i  +  1));  %  Determine  the  rate  needed  for  this  transition 

%ratesall ( runnum, i+1 )  =  rates(itl); 
i = i + 1  ; 

end 

%  Determine  the  number  of  transitions  from  one  state  to  the  next 

numij  =  zeros  (length (S) , length (S) ) ; 
for  k  =  1 : length  (x) -2 

numi j (x (k) , x (k+1 ) )  =  numi j (x (k) , x (k+1) )  +  1; 

end 

%  Compute  the  observed  transition  rates  among  the  environment  states 

O, 

for  i  =  1 : length (numi j ) 

for  j  =  1 : length (numi j ) 

if  j  ~=  i 

Qhat(i,j)  =  numi j ( i, j ) /sum ( statetimeTi j ( i ,:)) ; 

else 

Qhat ( i , j )  =  0  ; 

end 

end 
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Qhat  (i, i)  =  -sum (Qhat (i,  :)); 

statetimeTi ( i , 1 )  =  sum ( statetimeTi j (i, : ) ) ;  %  total  time  spent  in  state  i 

end 

o, 

%  This  section  creates  vectors  of  the  holding  times  in  States  1,2,  and  3 
%  ****IF  ADDITIONAL  STATES  ARE  ADDED,  THEY  MUST  BE  HARDCODED  HERE**** 

O, 

for  i  =  1 : length (time) 
for  j  =  1: length (S) 

if  rates  (i)  ==  R(j) 

HoldingTime ( i , j )  =  time(i);  %  creates  matrix  where  column  i 
end  %  represents  holding  times  in  state  i 

end 

end 

Stateltimes  =  HoldingTime (:, 1) ; 

Stateltimes  =  Stateltimes (find (Stateltimes) ) ;  %  eliminates  zeros  in  holding  time 

%  column 

State2times  =  HoldingTime (:, 2 ) ; 

State2times  =  State2times (find (State2times) ) ;  %  eliminates  zeros  in  holding  time 

%  column 

State3times  =  HoldingTime (:, 3) ; 

State3times  =  State3times (find (State3times) ) ;  %  eliminates  zeros  in  holding  time 

%  column 

o, 

%  This  section  computes  the  phase-type  approximations  for  each  of  the 
%  environment  states 

%  ****IF  ADDITIONAL  STATES  ARE  ADDED,  THEY  MUST  BE  HARDCODED  HERE**** 

O, 

M=10000;  %  this  is  the  large  value  for  the  'hot  potato  method' 

%  i.e.  since  there  is  really  no  'absorbing  phase'  for  any 
%  state,  we  wish  to  visit  each  absorbing  phase 
%  instantaneously. 

%  therefore,  its  rate  is  extremely  large  (approaches  infinity) 


for  State  =  1: length (S) 
clear  A; 
clear  alphal; 
clear  alpha2; 
clear  alpha; 
clear  beta; 
clear  k; 
clear  T; 
clear  To; 
if  State  ==  1 

distribution  =  2; 

A  =  Stateltimes; 

alphal=0;  %  parameter  not  needed;  dummy  value  only 
alpha2=0;  %  parameter  not  needed;  dummy  value  only 
alpha=K (State, 2 ) ; 
beta=K (State, 1 ) ; 

Aparam  =0;  %  parameter  not  needed;  dummy  value  only 
Bparam  =0;  %  parameter  not  needed;  dummy  value  only 
[k,  T, To]  = 

phasetypeapproxmixed (A, alpha, beta, alphal , alpha2 , Aparam, Bparam, distribution) ; 
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%  creates  phase-type  approximations  to  each  of  the  holding  time  distributions 
kl=k; 

T1=T; 

Tol=To; 

zer ovect or 1  =  zeros ( 1 ,  kl )  ; 

Psimatrixll  =  [  T1  Tol;  %  subgenerator  matrix  for  state 

zerovectorl  Qhat ( 1 , 1 ) *M] ; 
elseif  State  ==  2 

distribution  =  2; 

A  =  State2times; 

alphal=0;  %  parameter  not  needed;  dummy  value  only 
alpha2=0;  %  parameter  not  needed;  dummy  value  only 
alpha=K (State, 2) ; 
beta=K (State , 1 ) ; 

Aparam  =0;  %  parameter  not  needed;  dummy  value  only 
Bparam  =0;  %  parameter  not  needed;  dummy  value  only 
[k,T,To]  = 

phasetypeapproxmixed (A, alpha, beta, alpha 1 , alpha2 , Aparam, Bparam, distribution) ; 

%  creates  phase-type  approximations  to  each  of  the  holding  time  distributions 
k2=k; 

T2=T; 

To2=To; 

zer ovect or 2  =  zeros ( 1 ,  k2  )  ; 

Psimatrix22  =  [  T2  To2;  %  subgenerator  matrix  for  state  2 

zerovector2  Qhat (2 , 2 ) *M] ; 
elseif  State  ==  3 

distribution  =  2; 

A  =  State3times; 
alpha=K (State, 2) ; 
beta=K (State , 1 ) ; 

alphal=0;  %  parameter  not  needed;  dummy  value  only 
alpha2=0;  %  parameter  not  needed;  dummy  value  only 
Aparam  =0;  %  parameter  not  needed;  dummy  value  only 
Bparam  =0;  %  parameter  not  needed;  dummy  value  only 
[k,T,To]  = 

phasetypeapproxmixed (A, alpha, beta, alpha 1 , alpha2 , Aparam, Bparam, distribution) ; 

%  creates  phase-type  approximations  to  each  of  the  holding  time  distributions 
k3=k; 

T  3=T  ; 

To3=To; 

zer ovect or 3  =  zeros ( 1 ,  k3 )  ; 

Psimatrix33  =  [  T3  To3;  %  subgenerator  matrix  for  state  3 

zerovector3  Qhat ( 3 , 3 ) *M] ; 


end 


end 


1 


%  OFF-DIAGONAL  subgenerator  matrices 

%  ****IF  ADDITIONAL  STATES  ARE  ADDED,  THEY  MUST  BE  HARDCODED  HERE**** 

Psimatrixl2  =  zeros (kl+1 , k2+l) ; 

Psimat rix!2 (kl+1 , 1 ) =Qhat ( 1 , 2 ) *M; 


Psimatrixl3  =  zeros (kl+1 , k3+l) ; 
Psimat rix 13 (kl  +  1 , 1 ) =Qhat ( 1 , 3 ) *M; 

Psimatrix21  =  zeros (k2+l , kl+1); 


A-8 


Psimatrix21  (k.2  +  1 ,  1 )  =Qhat  (2 ,  1 )  *M; 

o, 

Psimatrix23  =  zeros  (k.2  +  1 ,  k3  +  l )  ; 
Psimatrix23 (k2+l , 1 ) =Qhat (2 , 3 ) *M; 

o, 

Psimatrix31  =  zeros (k3  +  l , kl  +  1 )  ; 
Psimatrix31 (k3+l , 1 ) =Qhat (3, 1 ) *M; 

o, 

Psimatrix32  =  zeros (k3  +  l , k2  +  l )  ; 
Psimatrix32 (k3+l , 1 ) =Qhat (3,2) *M; 


Psimatrix 


[Psimatrixll 
Psimatrix2 1 
Psimatrix31 


Psimatrixl2 

Psimatrix22 

Psimatrix32 


Psimatrixl3 ; 
Psimatrix23; 
Psimatrix33 ] ; 


%  overall  generator  matrix 
%  of  (newly-created)  CTMC 


Sr'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'te'k'k'k'k'k'k'k'k'k'k'k 


%  This  section  creates  the  expanded  degradation  rate  matrix  Lambda 

%  ****IF  ADDITIONAL  STATES  ARE  ADDED,  THEY  MUST  BE  HARDCODED  HERE**** 


V=diag (R) ; 

Lll  =  diag (diag (V (1, 1) *ones  (kl  +  1) )) ; 

L12  =  zeros (kl+1 , k2+l ) ; 

L13  =  zeros (kl+1, k3+l) ; 

L21  =  zeros (k2+l , kl+1 ) ; 

L22  =  diag (diag (V (2, 2) *ones  (k2  +  l) ) )  ; 

L23  =  zeros (k2+l , k3+l ) ; 

L31  =  zeros (k3+l, kl+1) ; 

L32  =  zeros (k3+l , k2+l ) ; 

L33  =  diag (diag (V (3, 3) *ones (k3+l) ) ) ; 

o, 

L  =  [Lll  L12  L13;  %  expanded  degradation  rate  matrix 

L21  L22  L23; 

L31  L32  L33] ; 

o, 

%  Display  the  results  of  the  phase-type  approximations 

o, 

fprintf  (  '  Size  of  generator  matrix  (Psimatrix):  %g  x  %g  \n',  length (Psimatrix) , 
length (Psimatrix) ) 
fprintf ( ' \n' ) 

fprintf  (  '  Size  of  expanded  degradation  rate  matrix  (Lambda)  :  %g  x  %g  \n', 
length  (L),  length  (L) ) 
fprintf ( ' \n ' ) 

fprintf ( '  ******************************  \n  '  ) 
fprintf ( ' \n ' ) 

o, 

%  This  section  computes  and  plots  the  failure  time  distribution  of  the 
%  system  exposed  to  the  SMP  environment 

o, 

for  w  =  l:length(t) 

FailureTimeProb (w)  =  invt_lap_Solo (t (w) , Psimatrix,  L) ; 
if  FailureTimeProb (w)  >  1 

FailureTimeProb (w)  =1;  %  eliminates  CDF  values  >  1 

elseif  FailureTimeProb (w)  <  0 

FailureTimeProb (w)  =0;  %  eliminates  CDF  values  <  0 

end 
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end 

o, 

^^^^^^iz^^iz^^iz^iz^^^^iz-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k'k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k'k-k-k-k-k-k-k-k-k-k-k-k-k 

%  Compute  first  and  second  moment  of  failure  time  distribution 

o, 

First_Moment  =  invt_lap_f irst_moment (WearThreshold, Psimatrix, L) ; 

o, 

Second_Moment  =  invt_lap_second_moment (WearThreshold, Psimatrix, L) ; 

o, 

Qriz^^iz^iz^^iz^^^^^^^iz^iziz-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k'k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

%  WRITE  OUTPUTS  TO  A  TAB  DELIMITED  TEXT  FILE 

O, 

fid  =  f open ( ' f : \AFIT\Thesis\Solo  Thesis\MATLAB  filesXSMP 
simulation\results\brakeanalytical . out ' , ' w ' ) ; 
for  k=l : length  (t ) 

fprintf (fid,  '%4.6f\t%4.6f\n',t (k) , FailureTimeProb (k) ) ; 

end 

f close (fid) ; 
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Appendix  B.  Code  for  Turbine  Blade  Example 


%  Plot_SMP_mixed.m 

%  This  MATLAB  program  runs  the  semi-Markov  process  (SMP)  simulation 
%  (semi_Markov.m)  and  then  the  'simulation  plus  phase-type  approximations  plus 
%  analytical  solution'  program  SMPmixedstates . m .  The  CDF  values  are  then 
%  compared  in  a  plot  (Figure  2) . 

%  Author:  Captain  Christopher  J.  Solo,  M.S.  candidate.  Operations  Research 
%  Air  Force  Institute  of  Technology 

%  Date:  January  29,  2004 


clear  all; 


t=0 . 005 : 0 . 05 : 50 . 0;  %  vector  of  time  points  at  which  to  evaluate  both  CDFs 

semi_Markov_mixed;  %  calls  program  to  simulate  semi-Markov  process  (SMP)  and 
obtain 

failure  times 


SMPmixedstates ; 


%  calls  program  to  simulate  SMP  (long  term),  compute  phase- 
%  type  approximations,  and  produce  analytical  results  of 
%  failure  time  distribution 


%  F  =  cdf  values  from  simulation 

%  FailureTimeProb  =  cdf  values  from  phase-type  analytical  solution 


%  Produces  plot  of  simulated  failure-time  distribution  vs.  analytical 
%  CDF  derived  via  conversion  of  environment  process  from  SMP  to  CTMC 

figure  (2)  ; 
plot (t, F, ' r ' ) ; 
hold  on; 

plot (t ,  FailureTimeProb,  ' b ' ) ; 
grid  off; 

title (' Failure-time  CDF  values:  Simulated  vs.  Analytical'); 
xlabel ( ' t ' ) ; 
ylabel ( ' P (T  <  t ) ' ) ; 

legend (' Simulated  CDF ',' Analytical  CDF',0); 
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%  PROGRAM  semi_Markov_mixed.m 

o, 

%  The  purpose  of  this  MATLAB  program  is  to  simulate  a  finite-state  semi-Markov 
%  process.  The  process  is  simulated  in  order  to  validate  analytical  solutions 
%  for  the  distribution  and  moments  of  the  lifetime  of  components  whose 
%  degradation  threshold  is  WearThreshold.  The  program  uses  function  "rando"  in 
%  order  to  select  the  next  state  after  a  state  transition. 

o, 

%  Orig  Author:  Jeffrey  P.  Kharoufeh,  Ph.D.  candidate,  IE  &  OR, 

%  Penn  State  University 

%  Date:  January  15,  2001 

%  Revised  by:  Captain  Chris  Solo,  M.S.  candidate,  OR,  \ 

%  Air  Force  Institute  of  Technology 

%  Date:  29  January  2004 

o, 

<^**************************************************************************** 

%  VARIABLE  DEFINTIONS 

*i(i(i(*i(ici(i(i(i(ici(icit-k-kic-k'k-k'k'k'ki('k'k'k-k'k'k-k'k-k'k-k'k'k'k-k'k-k'k-k'k'k'ki('k'k'k-k'k'k-k'k-k'k-k'k'k-k-kic-k 

o, 

%  WearThreshold  =  Degradation  threshold  (signifies  failure) 

%  k  =  an  index  variable 

%  VECTOR/FUNCTION  DEFINITIONS 

o, 

%  Lifetime  =  Vector  of  lifetimes 

%  B  =  Vector  for  the  initial  probability  distribution  of  the  Markov  process 
%  P  =  Probability  transition  matrix 

%  V  =  Vector  of  degradation  rates  (i.e.,  V(i)=  degradation  rate  when  environment 
%  is  in  state  i) . 

%  Z  =  Vector  of  environment  states 

%  The  program  assumes  a  state  space  of  the  form  S={1,2, . . .K} 

semi_Markov_input_mixed;  %  Obtain  intialization  parameters 

Z  =  [ ] ;  %  Initialize  Z  . 

for  k  =  1 : length (Lifetime) 

Z  =  El; 

Z ( 1 ) =rando (B) ;  %  Initial  state  of  the  environment  at  time  0 

if  (Z (1)  ==  S (1)  )  |  (Z (1)  ==  S (2)  ) 

Lifetime(k)  =  betarnd (K ( Z ( 1 ) , 1 ) , K ( Z ( 1 ) , 2 ) ) ;  %  Time  spent  in  initial  state 
elseif  (Z (1 )  -  S (3)  )  I  (Z ( 1)  ==  S  (4 )  ) 

Lifetime(k)  =  weibrnd (K ( Z ( 1 ) , 1 ) , K ( Z ( 1 ) , 2 ) ) ;  %  Time  spent  in  initial  state 

end 

D  =  0.0;  %  At  time  0,  degradation  is  0 

D  =  D  +  Lifetime (k) *V(Z (1) ) ;  %  Add  degradation  to  cumulative  degradation 

1  =  1; 

while  D  <  WearThreshold  %  Do  while  cumulative  degradation  is  <  WearThreshold 
Z(i+1)  =  rando (P (Z ( i ),:)) ;  %  Use  the  matrix  P  to  determine  next  state 

if  ( Z  ( i  + 1 )  ==  S  (1) )  |  (Z(i  +  1)  ==  S  (2 ) ) 

new_time  =  betarnd (K ( Z ( i  +  1 ), 1 ), K ( Z  ( i  +  1 ), 2 )) ;  %  Time  spent  in  state  i  +  1 
elseif  ( Z ( i+1 )  ==  S (3) )  I  (Z (i+1)  ==  S(4)) 
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new_time  =  weibrnd (K ( Z (i  +  1 ) , 1 ) , K ( Z  ( i  +  1 ) , 2 ) ) ;  %  Time  spent  in  state  i  +  1 

end 

D  =  D  +  new_time*V (Z (i+1) ) ;  %  Update  cumulative  degradation 
Lifetime (k)  =  Lifetime (k)  +  new_time;  %  Update  total  lifetime 

i  =  i  + 1  ; 

end 

Lifetime (k) =Lifetime (k) - (1/V (Z (i) ))* (D-WearThreshold) ;  %  Subtract  time  after 

%  reaching  threshold 

end 

get_cdf_mixed;  %  Computes  empirical  distribution 
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PROGRAM  semi_Markov_input_mixed . m 

The  purpose  of  this  MATLAB  program  is  to  provide  input  data  to  program 
semi_Markov_mixed . m  for  simulation  of  a  finite-state  semi-Markov  process  (K 
states) .  The  process  is  simulated  in  order  to  validate  analytical  solutions 
for  the  distribution  and  moments  of  the  failure  time  of  systems  with 
degradation  threshold  L  subject  to  a  SMP  environment.  The  program  uses 
function  "rando"  in  order  to  select  the  next  state  after  a  state  transition. 


Orig  Author: 

Date  : 
Revised  by: 

Last  Revised: 


Jeffrey  P.  Kharoufeh,  Ph.D.  candidate,  IE  &  OR, 
Penn  State  University 
January  15,  2001 

Captain  Christopher  Solo,  M.S.,  OR, 

Air  Force  Institute  of  Technology 
January  17,  2004 


VARIABLE  DEFINTIONS 


ifitic-k-k-k'k'k-k-k'k'k-k-k-k'k'k'k'k-k'k-k-k-k-k-k'k-k-k'k'k'k-k-k'k'k'k-k'k'k-k-k-k'k'k'k-k-k'k-k-k-k'k'k'k'k-k'k-k-k-k-k'k'k'k 


%  WearThreshold  =  degradation  threshold  (system  fails  after  accumulating  this 
%  amount  of  damage) 

o, 

%  VECTOR/FUNCTION  DEFINITIONS 

o, 

%  T  =  Vector  of  failure  time  observations 

%  B  =  Vector  for  the  initial  probability  distribution  of  the  semi-Markov  process 
%  P  =  Probability  transition  matrix 

%  V  =  Vector  of  degradation  rates  (i.e.,  V(i)=  degrad,  rate  when  environment  is 
%  in  state  i) . 

%  Initialize  variable  and  vector  values 


WearThreshold  =  20.0;  %  Component  reaches  'failure'  with  this  level  of 

%  degradation 

N  =  4;  %  number  of  states 

Lifetime  =  zeros  (1,  100000) ;  %  number  of  simulations  (failure  time  obs) 

B  =  [1  zeros (1, N-l )] ; 


s  = 

[1, 

.2,3,4] 

;  %  states  of  the  environment 

V  = 

o. 

[1. 

.7/10, 

4.3/10,  7.5/10,  12.9/10];  %  vector  of  degradation  rates 

p  = 

(0 

.6  .2 

.2;  %  transition  probability  matrix; 

.  5 

0  .4 

.  1; 

.2 

.4  0 

■  4; 

o. 

.  1 

.3  .6 

0]  ; 

K  = 

(3 

5; 

%  Parameter  matrix — each  row  gives  parameters  of  distribution 

2 

4; 

5 

2; 

6 

3]  ; 
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PROGRAM  get_cdf_mixed.m 

The  purpose  of  this  MATLAB  program  is  to  construct  an  empirical  distribution 
function  based  upon  a  vector  of  observations  of  some  continuous  random 
variable.  The  values  of  the  CDF  are  generally  used  for  the  purpose  of 
performing  hypothesis  test  on  the  equality  of  two  nonparametric  distributions. 
Other  approaches  for  testing  may  be  the  Cramer  von  Mises  test,  which  is 
similar  to  the  K-S  two-sample  test,  but  slightly  more  robust. 


Orig  Author: 

Date : 
Revised  by: 

Last  Revised: 


Jeffrey  P.  Kharoufeh,  Ph.D.  candidate,  IE  &  OR, 
Penn  State  University 
June  2000 

Captain  Chris  Solo,  M.S.  candidate,  OR, 

Air  Force  Institute  of  Technology 
January  18,  2004 


VECTOR  DEFINITIONS 


Lifetime  =  The  finite-dimensional  vector  of  real  obs 
t  =  The  vector  of  points  at  which  to  estimate  the  CDF. 

F  =  The  vector  of  CDF  estimates.  That  is,  F(i)  =  F(t(i)),  i=l , 2 , . . . , nof 

{1,2, . . ,,C} . 

★  ★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★it*********************** 


F  =  zeros ( 1 , length (t )) ;  %  This  ensures  that  t  and  F  are  vectors  of  equal  length 

O, 

%  Compute  the  cdf  value  at  the  point  tO 

o, 

for  i  =  1: length (t) 

F  ( i )  =  0.0; 

for  j  =  1 : length (Lifetime ) 
if  Lifetime ( j ) <=  t(i) 

F(i)  =  F(i)  +  1/length (Lifetime) ; 

else 

F  ( i )  =  F  ( i )  ; 

end 

end 

end 


WRITE  OUTPUTS  TO  A  TAB  DELIMITED  TEXT  FILE 


fid  =  fopen ( ' F : \AFIT\Thesis\Solo  Thesis\SMP 
simulation\results\turbinesimulated . out ' , ' w ' ) ; 
for  b=l : length (t ) 

fprintfffid,  '%4.6f\t%4.6f\n',t (b) ,F(b) ) ; 

end 

f close (fid) 
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%  SMPmixedstates . m 

O, 

%  The  purpose  of  this  MATLAB  program  is  to  simulate  a  semi-Markov  process  (SMP) 

%  on  a  state  space  S  and  convert  it  into  a  continuous-time  Markov  chain,  using  a 
%  phase-type  distribution  to  approximate  the  holding  time  distribution  of  each 
%  environment  state.  Each  of  the  environment  states  is  known  to  have  a  Weibull 
%  holding  time  distribution  (each  with  different  parameters.)  The  output 
%  is  the  analytical  solution  for  the  failure-time  distribution  of  a  system 
%  subject  to  the  SMP  environment,  based  on  the  results  of  Kharoufeh  (2003)  and 
%  Kharoufeh  and  Sipe  (2004)  . 

o, 

%  Original  code  for  process. m  and  rando.m  obtained  from  Craig  L.  Zirbel  at 
%  http://www-math.bgsu.edu/~zirbel/ap/  It  appeared  that  the  code  was  based 
%  on  a  discrete  time  Markov  chain,  and  Captain  Steve  Cox  modified  it  to  work 
%  for  a  continuous-time  Markov  chain.  Captain  Chris  Solo  modified  it  to  be 
%  a  SMP  (January  2004) . 

o, 

o, 

%  Author:  Captain  Christopher  J.  Solo,  M.S.  candidate, 

%  Operations  Research 

%  Air  Force  Institute  of  Technology 

%  Date:  September  13,  2003 

%  Last  Revised:  January  18,  2004 

o, 

%  USER  INPUTS 

%  Tmax  =  time  of  SMP  run 

%  numruns  =  number  of  simulations  (runs) 

%  mu  =  initial  probability  vector  of  environment  states 
%  S  =  state  space  vector 

%  R  =  degradation  rate  vector  of  environment  states 
%  P  =  transition  probability  matrix 

%  K  =  matrix  of  Weibull  parameters  for  environment  states 
%  M  =  large  value  used  in  Big  M  (hot  potato)  method 

o, 

%  Additionally,  rando.m,  phasetypeapproxmixed.m,  PHvalues.m,  invt_lap_Solo . m, 

%  and  failLST.m  are  needed. 

O, 

Sfe  ~k  -k  -k  "k  k  ★  k  k  k  ★  ★  k  ★  k  k  ★  k  k  k  k  k  ★  ★  ★  ★  k  k  k  k  ★  *  ★  k  k  k  k  ★  *  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  *  ~k  k  k  k  ★  *  *  k  *  k  k  k  ★  *  -k  k 

%  For  many  short  runs,  set  numruns  high  and  Tmax  low.  For  one  long  run, 

%  set  numruns  equal  to  1  and  make  Tmax  high. 

format  long  g; 

%  User  Inputs 

o, 

numruns  =  1;  %  the  number  of  times  (runs)  to  simulate  the  process 

Tmax  =  10000;  %  This  is  the  length  of  the  run  of  the  environment  process 

mu  =  [1,0,0,0];  %  Semi-Markov  process  (SMP)  starts  in  state  1 

S  =  [1,2,3,4];  %  There  are  three  states  defined 
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R  =  [1.7/10,  4.3/10,  7.5/10,  12.9/10];  %  Degradation  rate  of  each  state  of 

%  environment 

P  =  [0  .6  .2  .2;  %  transition  probability  matrix; 

.5  0.4.1; 

.2.4  0.4; 

.1.3.6  0]  ; 

K  =  [3  5;  %  Parameter  matrix — each  row  gives  parameters  of  distribution 

2  4; 

5  2; 

6  3]; 

O, 

O, 

%  Simulation  of  the  semi-Markov  process  (SMP) 

O, 

statetimeTij  =  zeros (length  (S) , length (S) ) ; 
statetimeTi  =  []; 
scale  =  1; 

O, 

TimeToFail ( 1 )  =0;  %  start  times  at  0 

%Tall ( runnum, 1 )  =  0; 

x (1)  =  rando (mu) ;  %  determine  starting  state  (time  0,  not  time  1) 

%xall (runnum, 1)  =  x(l); 
rates  (1)  =  R (x ( 1 ) ) ; 

%ratesall (runnum, 1)  =  rates  (1); 
i  =  1; 

o, 

while  TimeToFail ( i )  <  Tmax, 

x(i+l)  =  rando (P (x ( i ),:)) ;  %  Row  vector  from  P-matrix  determines  next 

%  transitions 

%xall  (runnum, i  +  1)  =  x(i  +  l); 
if  x (i)  ==  S ( 1 )  |  x(i)  ==  S (2 ) 

time(i)  =  betarnd (K(x(i) , 1) ,K(x(i) ,2) ) ;  %  Generate  beta  rv  based  on 

%  current  state 

elseif  x ( i )  ==  S(3)  I  x(i)  ==  S(4) 

time(i)  =  weibrnd (K (x ( i ) , 1 ) , K (x ( i ) , 2 ) ) ;  %  Generate  Weibull  rv  based  on 

%  current  state 

end 

statetimeTi j (x (i) , x (i+1) )  =  statetimeTi j (x ( i ), x ( i+1 ) )  +  time(i); 

TimeToFail (i+1 )  =  TimeToFail (i)  +  time(i);  %  Add  to  current  time 
%Tall (runnum, i+1)  =  T(i+1); 

rates  (i  +  1)  =  R(x(i  +  1));  %  Determine  the  rate  needed  for  this  transition 

%ratesall (runnum, i+1 )  =  rates (i+1); 
i=i+l ; 

end 

O, 

%  Determine  the  number  of  transitions  from  one  state  to  the  next 

O, 

numij  =  zeros (length (S) , length (S) ) ; 
for  k  =  1 : length (x) -2 

numi  j  (x  (k)  ,  x  (k+1)  )  =  numi  j  (x  (k)  ,  x  (k+1 )  )  +  1; 

end 

O, 

%  Compute  the  observed  transition  rates  among  the  environment  states 
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for  i  =  1 : length (numi j ) 

for  j  =  1 : length (numi j ) 
if  j  ~=  i 

Qhat(i,j)  =  numi j (i, j ) /sum (statetimeTi j (i,  :))  ; 

else 

Qhat  (i,  j)  =  0; 

end 

end 

Qhat  (i, i)  =  -sum (Qhat (i,  :)); 

statetimeTi ( i , 1 )  =  sum ( statetimeTi j (i, :)) ;  %  total  time  spent  in  state  i 

end 

o, 

%  This  section  creates  vectors  of  the  holding  times  in  States  1,2,  and  3 
%  ****IF  ADDITIONAL  STATES  ARE  ADDED,  THEY  MUST  BE  HARDCODED  HERE**** 

O, 

for  i  =  1 : length (time) 
for  j  =  1: length (S) 

if  rates  (i)  ==  R(j) 

HoldingTime ( i , j )  =  time(i);  %  creates  matrix  where  column  i 
end  %  represents  holding  times  in  state  i 

end 

end 

Stateltimes  =  HoldingTime (:, 1) ; 

Stateltimes  =  Stateltimes (find (Stateltimes) ) ;  %  eliminates  zeros  in  holding  time 

%  column 

State2times  =  HoldingTime (:, 2 ) ; 

State2times  =  State2times (find (State2times) ) ;  %  eliminates  zeros  in  holding  time 

%  column 

State3times  =  HoldingTime (:, 3) ; 

State3times  =  State3times (find (State3times) ) ;  %  eliminates  zeros  in  holding  time 

%  column 

State4times  =  HoldingTime (:, 4 ) ; 

State4times  =  State4times (find (State4times) ) ;  %  eliminates  zeros  in  holding  time 

%  column 

O, 

%  This  section  computes  the  phase-type  approximations  for  each  of  the 
%  environment  states 

%  ****IF  ADDITIONAL  STATES  ARE  ADDED,  THEY  MUST  BE  HARDCODED  HERE**** 

O, 

M=10000;  %  this  is  the  large  value  for  the  'hot  potato  method' 

%  i.e.  since  there  is  really  no  'absorbing  phase'  for  any 
%  state,  we  wish  to  visit  each  absorbing  phase 
%  instantaneously. 

%  therefore,  its  rate  is  extremely  large  (approaches  infinity) 


for  State  =  1: length (S) 
clear  A; 
clear  alphal; 
clear  alpha2; 
clear  alpha; 
clear  beta; 
clear  k; 
clear  T; 
clear  To; 
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if  State  ==  1 

distribution  =  3; 

A  =  Stateltimes; 
alphal=K (State, 1) ; 
alpha2=K (State, 2 ) ; 

alpha=0;  %  parameter  not  needed;  dummy  value  only 
beta=0;  %  parameter  not  needed;  dummy  value  only 
Aparam=0;  %  parameter  not  needed;  dummy  value  only 
Bparam=0;  %  parameter  not  needed;  dummy  value  only 
[k,  T,  To]  = 

phasetypeapproxmixed (A, alpha, beta, alphal , alpha2 , Aparam, Bparam, distribution) ; 

%  creates  phase-type  approximations  to  each  of  the  holding  time  distributions 
kl=k; 

Tl-'T; 

Tol=To; 

zerovectorl  =  zeros (1,  kl)  ; 

Psimatrixll  =  [  T1  Tol;  %  subgenerator  matrix  for  state  1 

zerovectorl  Qhat ( 1 , 1 ) *M] ; 

elseif  State  ==  2 

distribution  =  3; 

A  =  State2times; 
alphal=K (State, 1) ; 
alpha2=K (State,  2 )  ; 

alpha=0;  %  parameter  not  needed;  dummy  value  only 
beta=0;  %  parameter  not  needed;  dummy  value  only 
Aparam=0;  %  parameter  not  needed;  dummy  value  only 
Bparam=0;  %  parameter  not  needed;  dummy  value  only 
[  k , T , To ]  = 

phasetypeapproxmixed (A, alpha, beta, alphal , alpha2 , Aparam, Bparam,  distribution)  ; 

%  creates  phase-type  approximations  to  each  of  the  holding  time  distributions 
k2=k; 

T2=T; 

To2=To; 

zerovector2  =  zeros (1,  k2)  ; 

Psimatrix22  =  [  T2  To2;  %  subgenerator  matrix  for  state  2 

zerovector2  Qhat (2 , 2 ) *M] ; 
elseif  State  ==  3 

distribution  =  2; 

A  =  State3times; 
alpha=K (State, 2) ; 
beta=K (State, 1) ; 

alphal=0;  %  parameter  not  needed;  dummy  value  only 
alpha2=0;  %  parameter  not  needed;  dummy  value  only 
Aparam=0;  %  parameter  not  needed;  dummy  value  only 
Bparam=0;  %  parameter  not  needed;  dummy  value  only 
[  k , T , To ]  = 

phasetypeapproxmixed (A, alpha, beta, alphal , alpha2 , Aparam, Bparam,  distribution)  ; 

%  creates  phase-type  approximations  to  each  of  the  holding  time  distributions 
k3=k; 

T3=T; 

To3=To; 

zerovector3  =  zeros (1,  k3)  ; 

Psimatrix33  =  [  T3  To3;  %  subgenerator  matrix  for  state  3 

zerovector3  Qhat  ( 3 , 3 ) *M] ; 
elseif  State  ==  4 

distribution  =  2; 

A  =  State4times; 
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alpha=K (State, 2) ; 
beta=K (State , 1 ) ; 

alphal=0;  %  parameter  not  needed;  dummy  value  only 
alpha2=0;  %  parameter  not  needed;  dummy  value  only 
Aparam=0;  %  parameter  not  needed;  dummy  value  only 
Bparam=0;  %  parameter  not  needed;  dummy  value  only 
[k,T,To]  = 

phasetypeapproxmixed (A, alpha, beta, alphal , alpha2 , Aparam, Bparam, distribution) ; 

%  creates  phase-type  approximations  to  each  of  the  holding  time  distributions 
k4=k; 

T4=T  ; 

To4=To; 

zerovector4  =  zeros ( 1 ,  k4 )  ; 

Psimatrix44  =  [  T4  To4;  %  subgenerator  matrix  for  state 

zerovector4  Qhat ( 4 , 4 ) *M] ; 


end 


end 


3 


%  OFF-DIAGONAL  subgenerator  matrices 

%  ****IF  ADDITIONAL  STATES  ARE  ADDED,  THEY  MUST  BE  HARDCODED  HERE**** 

Psimatrixl2  =  zeros (kl+1 , k2+l) ; 

Psimat rixl2 (kl+1 , 1 ) =Qhat ( 1 , 2 ) *M; 

Psimatrixl3  =  zeros (kl+1 , k3+l) ; 

Psimat rix 13 (kl+1 , 1 ) =Qhat ( 1 , 3 ) *M; 

Psimatrixl4  =  zeros (kl+1 , k4+l) ; 

Psimat rix 14 (kl+1 , 1 ) =Qhat ( 1 , 4 ) *M; 

Psimatrix21  =  zeros (k2+l , kl+1) ; 

Psimat rix2 1 (k2+l , 1 ) =Qhat (2 , 1 ) *M; 

Psimatrix23  =  zeros (k2+l, k3+l) ; 

Psimat rix2 3 (k2  +  l , 1 ) =Qhat (2 , 3 ) *M; 

Psimatrix24  =  zeros (k2+l, k4+l) ; 

Psimat rix2 4 (k2+l , 1 ) =Qhat (2 , 4 ) *M; 

Psimatrix31  =  zeros (k3+l , kl+1) ; 

Psimatrix31 (k3  +  l , 1 ) =Qhat (3,  1 ) *M; 

Psimatrix32  =  zeros (k3+l , k2+l) ; 

Psimatrix32 (k3+l , 1 ) =Qhat (3, 2 ) *M; 

Psimatrix34  =  zeros (k3+l , k4+l) ; 

Psimatrix34 (k3  +  l , 1 ) =Qhat (3,  4 ) *M; 

Psimatrix41  =  zeros (k4+l , kl+1) ; 

Psimat rix 4 1 (k4+l , 1 ) =Qhat ( 4 , 1 ) *M; 

Psimatrix42  =  zeros (k4  +  l ,  k2  +  l)  ; 

Psimat rix 4 2 (k4+l , 1 ) =Qhat ( 4 , 2 ) *M; 

Psimatrix43  =  zeros (k4+l , k3+l) ; 

Psimat rix 4 3 (k4+l , 1 ) =Qhat ( 4 , 3 ) *M; 
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Psimatrix 


[Psimatrixll 
Psimatrix2 1 
Psimatrix31 
Psimatrix4 1 


Psimatrixl2 

Psimatrix22 

Psimatrix32 

Psimatrix42 


Psimatrixl3 

Psimatrix23 

Psimatrix33 

Psimatrix43 


Psimatrixl4 ; 
Psimatrix24 ; 
Psimatrix34 ; 
Psimatrix4  4 ] ; 


%  overall 
%  gener.  matrix 
%  of  (newly- 
%  created)  CTMC 


%  This  section  creates  the  expanded  degradation  rate  matrix  Lambda 

%  ****IF  ADDITIONAL  STATES  ARE  ADDED,  THEY  MUST  BE  HARDCODED  HERE**** 


V=diag (R) ; 

Lll  =  diag (diag (V  (1, 1) 

*ones (kl+1) ) ) ; 

L12 

=  zeros (kl+1 , k2+l ) 

r 

L13 

=  zeros  (kl  +  1 ,  k.3  +  1 ) 

r 

L14 

=  zeros  (kl  +  1 ,  k.4  +  1 ) 

r 

L21 

=  zeros  (k.2  +  1 ,  kl  +  1 ) 

r 

L22 

=  diag (diag (V (2, 2) 

*ones  (k2  +  l ) ) ) ; 

L23 

=  zeros  (k.2  +  1 ,  k.3  +  1 ) 

} 

L24 

=  zeros  (k.2  +  1 ,  k4  +  l ) 

r 

L31 

=  zeros (k3+l , kl+1 ) 

r 

L32 

=  zeros  (k.3  +  1 ,  k.2  +  1 ) 

r 

L33 

=  diag (diag (V (3,  3) 

*ones  (k3  +  l ) ) ) ; 

L34 

=  zeros  (k.3  +  1 ,  k4  +  l ) 

} 

L41 

=  zeros  (k.4  +  1 ,  kl  +  1 ) 

r 

L42 

=  zeros  (k.4  +  1 ,  k.2  +  1 ) 

r 

L43 

=  zeros  (k.4  +  1 ,  k3  +  l ) 

} 

L44 

o. 

=  diag (diag (V ( 4 , 4 ) 

*ones  (k4  +  l ) ) ) ; 

L  = 

[Lll  L12  L13  L14 ; 

%  expanded  degradation  rate  matrix 

L21  L22  L23  L24; 

L31  L32  L33  L34 ; 

L41  L42  L43  L44] ; 

O, 

%  Display  the  results 

of  the  phase-type  approximations 

fprintf ( '  Size  of  generator  matrix  (Psimatrix):  %g  x  %g  \n',  length (Psimatrix) , 
length (Psimatrix) ) 
fprintf ( ' \n ' ) 

fprintf ( '  Size  of  expanded  degradation  rate  matrix  (Lambda) :  %g  x  %g  \n', 
length  (L),  length  (L) ) 
fprintf ( ' \n ' ) 

fprintf ( '  ******************************  \n  '  ) 
fprintf ( ' \n' ) 

o, 

%  This  section  computes  and  plots  the  failure  time  distribution  of  the 
%  system  exposed  to  the  SMP  environment 

o, 

for  w  =  l:length(t) 

FailureTimeProb (w)  =  invt_lap_Solo (t (w) , Psimatrix, L) ; 
if  FailureTimeProb (w)  >  1 


FailureTimeProb (w) 
elseif  FailureTimeProb 

=  1; 
(w)  < 

%  eliminates 

0 

CDF 

values 

>  1 

FailureTimeProb (w) 

end 

=  0; 

%  eliminates 

CDF 

values 

<  0 

end 
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%  Compute  first  and  second  moment  of  failure  time  distribution 

o, 

First_Moment  =  invt_lap_f irst_moment (WearThreshold, Psimatrix, L) ; 

o, 

Second_Moment  =  invt_lap_second_moment (WearThreshold, Psimatrix, L) ; 

o, 

2r^^^^^iz^^iz^^iz^-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k'k-k-k-k-k-k-k-k-k-k-k-k-k 

o, 

%  WRITE  OUTPUTS  TO  A  TAB  DELIMITED  TEXT  FILE 

O, 

fid  =  f open ( ' F : \AFIT\Thesis\Solo  ThesisXSMP 
simulation\results\turbineanalytical . out ' , ' w ' ) ; 
for  k=l : length  (t) 

fprintf (fid,  '%4.6f\t%4.6f\n',t (k) , FailureTimeProb (k) ) ; 

end 

f close (fid) ; 
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Appendix  C.  Code  for  Coating  Example 


%  Plot_SMP_f ive .m 

o, 

%  This  MATLAB  program  runs  the  semi-Markov  process  (SMP)  simulation 
%  ( semi_Markov_f ive . m)  and  then  the  'simulation  plus  phase-type  approximations 
%  plus  analytical  solution'  program  SMPfivestates .m.  The  CDF  values  are  then 
%  compared  in  a  plot  (Figure  2) . 

o, 

O, 

%  Author:  Captain  Christopher  J.  Solo,  M.S.  candidate,  Operations  Research 
%  Air  Force  Institute  of  Technology 

%  Date:  January  29,  2004 

o, 

o, 

clear  all; 


t=0 .05:0 .005 : 10; 

o, 

semi_Markov_f ive; 
obtain 

o, 

o, 

SMPfivestates ; 


%  vector  of  time  points  at  which  to  evaluate  both  CDFs 

%  calls  program  to  simulate  semi-Markov  process  (SMP)  and 
failure  times 

%  calls  program  to  simulate  SMP  (long  term),  compute  phase-type 
%  approximations,  &  produce  analytical  results  of  failure  time 
%  distribution 


%  F  =  cdf  values  from  simulation 

%  FailureTimeProb  =  cdf  values  from  phase-type  analytical  solution 

o, 

%  Produces  plot  of  simulated  failure-time  distribution  vs.  analytical 
%  CDF  derived  via  conversion  of  environment  process  from  SMP  to  CTMC 


figure  (2) ; 
plot (t , F, ' r ' ) ; 
hold  on; 

plot (t ,  FailureTimeProb,  ' b ' ) ; 
grid  off; 

title (' Failure-time  CDF  values:  Simulated  vs.  Analytical'); 
xlabel ( ' t '  )  ; 
ylabel ( 'P (T  <  t)  '  )  ; 

legend (' Simulated  CDF ',' Analytical  CDF',0); 
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PROGRAM  semi_Markov_f ive . m 

The  purpose  of  this  MATLAB  program  is  to  simulate  a  finite-state  semi-Markov 
process.  The  process  is  simulated  in  order  to  validate  analytical  solutions 
for  the  distribution  and  moments  of  the  failure  time  of  systems  with 
degradation  threshold  L  subject  to  a  SMP  environment.  The  program  uses 
function  "rando"  in  order  to  select  the  next  state  after  a  state  transition. 


Orig  Author : 

Date : 
Revised  by: 

Last  Revised: 

■k'k'k'k-k'k'k-k'k'k'k-k'k-k 


Jeffrey  P.  Kharoufeh,  Ph.D.  candidate,  IE  &  OR, 
Penn  State  University 
January  15,  2001 

Captain  Christopher  Solo,  M.S.,  OR, 

Air  Force  Institute  of  Technology 
January  18,  2004 


VARIABLE  DEFINTIONS 


WearThreshold  =  degradation  threshold  (system  fails  after  accumulating  this 
amount  of  damage) 
k  =  an  index  variable 

★★★★★it************************************************************************* 


%  VECTOR/FUNCTION  DEFINITIONS 

O, 

%  T  =  Vector  of  failure  time  observations 

%  B  =  Vector  for  the  initial  probability  distribution  of  the  semi-Markov  process 
%  P  =  Probability  transition  matrix 

%  V  =  Vector  of  degradation  rates  (i.e.,  V(i)=  degrad,  rate  when  environment  is 
%  in  state  i) . 

%  Z  =  Vector  of  environment  states 

o, 

%  The  program  assumes  a  state  space  of  the  form  S={1,2, . . .K} 

o, 

%  Simulation  of  semi-Markov  process  (SMP) 


semi_Markov_input_f ive;  %  Obtain  initialization  parameters 

Z  =  [];  %  Initialize  Z. 


for  k  =  1 : length (Lifetime) 

Z  =  El; 

Z ( 1 ) =rando (B) ; 


if  Z  ( 1 )  ==  S  (1 )  I  Z  ( 1 )  ==  S  (3) 

Lifetime(k)  =  betarnd (K ( Z ( 1 ) , 1 ) , K ( Z ( 1 ) , 2 ) ) ;  %  Time  spent  in  initial  state 
elseif  Z ( 1 )  ==  S (2 )  I  Z  (1 )  ==  S  (4 ) 

Lifetime(k)  =  weibrnd (K ( Z ( 1 ) , 1 ) , K ( Z ( 1 ) , 2 ) )  ; 
elseif  Z (1)  ==  S (5) 

Lifetime(k)  =  gamrnd (K ( Z ( 1 ) , 1 ) , K ( Z ( 1 ) , 2 ) ) ; 

end 

D  =  0.0;  %  At  time  0,  distance  traveled  is  0 

D  =  D  +  Lifetime (k) *V(Z (1) ) ;  %  Add  traveled  distance  to  cumulative  distance 

i  =  l; 

while  D  <  WearThreshold  %  Do  while  distance  traveled  is  less  than  WearThr 
Z(i+1)  =  rando (P (Z ( i ),:)) ;  %  Use  the  matrix  P  to  determine  next  state 

if  Z  (i+1)  ==  S ( 1 )  I  Z  (i  +  1)  ==  S  (3) 
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new_time  =  betarnd (K (Z (i  +  1) , 1) , K (Z (i  +  1 ) , 2 ) )  ; 

%  Generate  beta  rv  based  on  current  state 

elseif  Z  (i+1)  ==  S(2)  |  Z (i  +  1)  ==  S(4) 

new_time  =  weibrnd (K (Z (i  +  1) , 1) , K  (Z (i  +  1 ) , 2 ) ) ; 

%  Generate  beta  rv  based  on  current  state 
elseif  Z (i+1)  ==  S  (5 ) 

new_time  =  gamrnd (K ( Z ( i+1 ) , 1 ) , K ( Z (i+1 ) , 2 ) ) ; 

%  Generate  gamma  rv  based  on  current  state 
end 

D  =  D  +  new_time*V (Z (i+1) ) ;  %  Update  cumulative  damage 

Lifetime (k)  =  Lifetime (k)  +  new_time;  %  Update  total  lifetime 
i=i  +  l 

end 

Lif et ime (k) =Lifetime (k) - ( 1/V ( Z ( i ) ) ) * (D-WearThreshold) ;  %  Subtract  extra  time 

%  after  WearThr  is 
%  reached 

end 

get_cdf_f ive;  %  computes  the  empirical  distribution  of  observed  failure  times 
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PROGRAM  semi_Markov_input_f ive . m 

The  purpose  of  this  MATLAB  program  is  to  provide  input  data  to  program 
semi_Markov_f ive . m  for  simulation  of  a  finite-state  semi-Markov  process  (K 
states) .  The  process  is  simulated  in  order  to  validate  analytical  solutions 
for  the  distribution  and  moments  of  the  failure  time  of  systems  with 
degradation  threshold  L  subject  to  a  SMP  environment.  The  program  uses 
function  "rando"  in  order  to  select  the  next  state  after  a  state  transition. 


Orig  Author: 

Date  : 
Revised  by: 

Last  Revised: 


Jeffrey  P.  Kharoufeh,  Ph.D.  candidate,  IE  &  OR, 
Penn  State  University 
January  15,  2001 

Captain  Christopher  Solo,  M.S.,  OR, 

Air  Force  Institute  of  Technology 
January  17,  2004 


iticiticititificitititicltitititicicjtiticltitititicicjtiticicieititicicjtititiclt'k-k'k-k-k-kit'k-k-k-k-k-k-k-kit'k-k-k-k-k-k-k-kit-k'k-k'k-k-k-k-k-k-k'k-k'k 

VARIABLE  DEFINTIONS 


%  WearThreshold  =  degradation  threshold  (system  fails  after  accumulating  this 
%  amount  of  damage) 

o, 

%  VECTOR/FUNCTION  DEFINITIONS 

o, 

%  T  =  Vector  of  failure  time  observations 

%  B  =  Vector  for  the  initial  probability  distribution  of  the  semi-Markov  process 
%  P  =  Probability  transition  matrix 

%  V  =  Vector  of  degradation  rates  (i.e.,  V(i)=  degrad,  rate  when  environment  is 
%  in  state  i) . 

%  Initialize  variable  and  vector  values 


WearThreshold  =  5.0; 

N  =  5;  %  number  of  states 

Lifetime  =  zeros  (1,  100000) ;  %  number  of  simulations  (failure  time  obs) 

B  =  [1  zeros (1, N-l )] ; 


s  = 

[1 

2  3 

4 

5]  ; 

%  states  of  the  environment 

V  = 

o. 

(0. 

.  4  6 

0. 

82  1 

.10  1.34  1.98];  %  matrix  of  degradation  rates 

p  = 

(0 

.5 

.2 

.2 

.1;  %  transition  probability  matrix; 

.  5 

0 

.3 

.  1 

.1; 

.2 

.  4 

0 

.2 

.2; 

.  1 

.  1 

.2 

0 

.6; 

.2 

.  1 

.  1 

.  6 

0]  ; 

K  = 

(3 

5; 

O, 

Parameter  matrix — each  row  gives  parameters  of  distribution 

5 

2; 

4 

6; 

6 

3 

.  5 

.10 

] ; 
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★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★ 
PROGRAM  get_cdf_f ive . m 

The  purpose  of  this  MATLAB  program  is  to  construct  an  empirical  distribution 
function  based  upon  a  vector  of  observations  of  some  continuous  random 
variable.  The  values  of  the  CDF  are  generally  used  for  the  purpose  of 
performing  hypothesis  test  on  the  equality  of  two  nonparametric  distributions. 
Other  approaches  for  testing  may  be  the  Cramer  von  Mises  test,  which  is 
similar  to  the  K-S  two-sample  test,  but  slightly  more  robust. 


Orig  Author: 

Date : 
Revised  by: 

Last  Revised: 


Jeffrey  P.  Kharoufeh,  Ph.D.  candidate,  IE  &  OR, 

Penn  State  University 
June  2000 

Captain  Chris  Solo,  M.S.  candidate,  OR, 

Air  Force  Institute  of  Technology 
January  18,  2004 

★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★it** 


VECTOR  DEFINITIONS 

★  ★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★it******** 


Lifetime  =  The  finite-dimensional  vector  of  real  observations 
t  =  The  vector  of  points  at  which  to  estimate  the  CDF. 

F  =  The  vector  of  CDF  estimates.  That  is,  F(i)  =  F(t(i)),  i  =  l ,  2 ,  . . . , nof 

{1,2,  .  .  ,,C}  . 

★  ★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★it**** 

=  zeros ( 1 , length (t ))  ;  %  This  ensures  that  t  and  F  are  vectors  of  equal  length 

Compute  the  cdf  value  at  the  point  tO 


for  i  =  1: length (t) 

F  ( i )  =  0.0; 

for  j  =  1 : length (Lifetime ) 
if  Lifetime ( j ) <=  t(i) 

F(i)  =  F(i)  +  1/length (Lifetime) ; 

else 

F  (i)  =  F  ( i )  ; 

end 

end 

end 

O, 

%  WRITE  OUTPUTS  TO  A  TAB  DELIMITED  TEXT  FILE 

O, 

fid  =  fopen ( ' F : \AFIT\Thesis\Solo  Thesis\SMP 
simulation\results\chemsimulated . out ' , ' w ' ) ; 
for  b=l : length (t ) 

fprintf(fid,  '%4.6f\t%4.6f\n',t (b) ,F(b) ) ; 

end 

f close (fid) 
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%  SMPf ivestates  .m 

o, 

%  The  purpose  of  this  MATLAB  program  is  to  simulate  a  semi-Markov  process  (SMP) 

%  on  a  state  space  S  and  convert  it  into  a  continuous-time  Markov  chain,  using  a 
%  phase-type  distribution  to  approximate  the  holding  time  distribution  of  each 
%  environment  state.  Each  of  the  environment  states  is  known  to  have  a  Weibull 
%  holding  time  distribution  (each  with  different  parameters.)  The  output 
%  is  the  analytical  solution  for  the  failure-time  distribution  of  a  system 
%  subject  to  the  SMP  environment,  based  on  the  results  of  Kharoufeh  (2003)  and 
%  Kharoufeh  and  Sipe  (2004)  . 

o, 

%  Original  code  for  process. m  and  rando.m  obtained  from  Craig  L.  Zirbel  at 
%  http://www-math.bgsu.edu/~zirbel/ap/  It  appeared  that  the  code  was  based 
%  on  a  discrete  time  Markov  chain,  and  Captain  Steve  Cox  modified  it  to  work 
%  for  a  continuous-time  Markov  chain.  Captain  Chris  Solo  modified  it  to  be 
%  a  SMP  (January  2004) . 

o, 

o, 

%  Author:  Captain  Christopher  J.  Solo,  M.S.  candidate, 

%  Operations  Research 

%  Air  Force  Institute  of  Technology 

%  Date:  September  13,  2003 

%  Last  Revised:  January  18,  2004 

o, 

%  USER  INPUTS 

%  Tmax  =  time  of  SMP  run 

%  numruns  =  number  of  simulations  (runs) 

%  mu  =  initial  probability  vector  of  environment  states 
%  S  =  state  space  vector 

%  R  =  degradation  rate  vector  of  environment  states 
%  P  =  transition  probability  matrix 

%  K  =  matrix  of  Weibull  parameters  for  environment  states 
%  M  =  large  value  used  in  Big  M  (hot  potato)  method 

o, 

%  Additionally,  rando.m,  phasetypeapproxmixed.m,  PHvalues.m,  invt_lap_Solo . m, 

%  and  failLST.m  are  needed. 

O, 

Sfe  ~k  -k  -k  "k  k  ★  k  k  k  ★  ★  k  ★  k  k  ★  k  k  k  k  k  ★  ★  ★  ★  k  k  k  k  ★  *  ★  k  k  k  k  ★  *  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  *  ~k  k  k  k  ★  *  *  k  *  k  k  k  ★  *  -k  k 

%  For  many  short  runs,  set  numruns  high  and  Tmax  low.  For  one  long  run, 

%  set  numruns  equal  to  1  and  make  Tmax  high. 

format  long  g; 

%  User  Inputs 

o, 

numruns  =  1;  %  the  number  of  times  (runs)  to  simulate  the  process 

Tmax  =  40000;  %  This  is  the  length  of  the  run  of  the  environment  process 

mu  =  [1,0, 0,0,0];  %  Semi -Markov  process  (SMP)  starts  in  state  1 

S  =  [1,2, 3, 4, 5];  %  There  are  three  states  defined 
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6.7,  9.9]; 


Degradation  rate  of  each  state  of  environment 


R  =  [2 

.3,  4. 

.1, 

5.5 

P  =  [0 

.5  .2 

.2 

.1; 

.  5 

0  .3 

.  1 

.1; 

.2 

.4  0 

.2 

.2; 

.  1 

.  1  .2 

0 

.6; 

.2 

.1  .1 

.  6 

0] 

K  =  [3 

5; 

distribution 

5 

2; 

4 

6; 

6 

3 

.  5 

q. 

-10]  ; 

%  Simulation 

of 

the 

transition  probability  matrix; 


Parameter  matrix — each  row  gives  parameters  of  different 


statetimeTij  =  zeros (length (S) , length (S) )  ; 
statetimeTi  =  []; 
scale  =  1; 

o, 

TimeToFail (1)  =0;  %  start  times  at  0 

%Tall (runnum, 1)  =  0; 

x(l)  =  rando (mu) ;  %  determine  starting  state  (time  0,  not  time  1) 

%xall ( runnum, 1 )  =  x(l); 
rates  (1)  =  R (x ( 1 ) ) ; 

%ratesall ( runnum,  1 )  =  rates  (1); 

i  =  1; 


while  TimeToFail (i)  <  Tmax, 

x(i+l)  =  rando (P (x ( i ),:)) ;  %  Row  vector  from  P-matrix 

determines  next  transitions 

%xall ( runnum, i+1 )  =  x(i+l); 
if  x  (i)  ==  S  (1 )  I  x  (i)  ==  S  (3) 

time(i)  =  betarnd (K (x (i) , 1) , K (x (i) , 2) )  ; 

%  Generate  beta  rv  based  on  current  state 
elseif  x ( i )  ==  S(2)  |  x(i)  ==  S  ( 4 ) 

time(i)  =  weibrnd (K (x (i) , 1) , K (x (i) , 2) ) ; 

%  Generate  Weibull  rv  based  on  current  state 
elseif  x (i)  ==  S ( 5 ) 

time(i)  =  gamrnd (K (x ( i ) , 1 ) , K (x (i )  ,  2 ) ) ; 

%  Generate  gamma  rv  based  on  current  state 
end 

statetimeTi j (x ( i ), x ( i+1 ) )  =  statetimeTi j (x ( i ), x ( i+1 ) )  +  time(i); 

TimeToFail ( i+1 )  =  TimeToFail (i)  +  time(i);  %  Add  to  current  time 
%Tall (runnum, i  +  1)  =  T ( i  +  1 )  ; 

rates  (i  +  1)  =  R(x(i  +  1));  %  Determine  the  rate  needed  for  this  transition 

%ratesall ( runnum, i+1 )  =  rates (i+1); 
i=i+l ; 

end 

o, 

%  Determine  the  number  of  transitions  from  one  state  to  the  next 


numij  =  zeros (length (S) , length (S) ) ; 
for  k  =  1 : length (x) -2 
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numi j (x (k) , x (k+1) )  =  numi j (x (k) , x (k+1 ) )  +  1; 

end 

O, 

%  Compute  the  observed  transition  rates  among  the  environment  states 

O, 

for  i  =  1 : length (numi j ) 

for  j  =  1 : length (numi j ) 
if  j  ~=  i 

Qhat(i,j)  =  numi j (i, j ) /sum (statetimeTi j (i,  :))  ; 

else 

Qhat  (i, j )  =  0 ; 

end 

end 

Qhat(i,i)  =  -sum (Qhat (i, :)) ; 

statetimeTi ( i , 1 )  =  sum ( statetimeTi j (i, :)) ;  %  total  time  spent  in  state  i 

end 

o, 

%  This  section  creates  vectors  of  the  holding  times  in  States  1,2,  and  3 
%  ****IF  ADDITIONAL  STATES  ARE  ADDED,  THEY  MUST  BE  HARDCODED  HERE**** 

O, 

for  i  =  1 : length (time) 
for  j  =  1: length (S) 

if  rates  (i)  ==  R(j) 

HoldingTime ( i , j )  =  time(i);  %  creates  matrix  where  column  i 
end  %  represents  holding  times  in  state  i 

end 

end 

Stateltimes  =  HoldingTime (:, 1) ; 

Stateltimes  =  Stateltimes (find (Stateltimes) ) ; 

%  eliminates  zeros  in  holding  time  column 
State2times  =  HoldingTime (:,  2 )  ; 

State2times  =  State2times (find (State2times) ) ; 

%  eliminates  zeros  in  holding  time  column 
State3times  =  HoldingTime (:, 3) ; 

State3times  =  State3times (find (State3times) ) ; 

%  eliminates  zeros  in  holding  time  column 
State4times  =  HoldingTime (:,  4 )  ; 

State4times  =  State4times (find (State4times) ) ; 

%  eliminates  zeros  in  holding  time  column 
State5times  =  HoldingTime (:, 5) ; 

State5times  =  State5times (find (State5times) ) ; 

%  eliminates  zeros  in  holding  time  column 

o, 

%  This  section  computes  the  phase-type  approximations  for  each  of  the 
environment  states 

%  ****IF  ADDITIONAL  STATES  ARE  ADDED,  THEY  MUST  BE  HARDCODED  HERE**** 

O, 

M=10000;  %  this  is  the  large  value  for  the  'hot  potato  method' 

%  i.e.  since  there  is  really  no  'absorbing  phase'  for  any 
%  state,  we  wish  to  visit  each  absorbing  phase 
%  instantaneously. 

%  therefore,  its  rate  is  extremely  large  (approaches  infinity) 

o, 

for  State  =  1: length (S) 
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clear  A; 
clear  alphal; 
clear  alpha2; 
clear  alpha; 
clear  beta; 
clear  Aparam; 
clear  Bparam; 
clear  k; 
clear  T; 
clear  To; 
if  State  ==  1 

distribution  =  3; 

A  =  Stateltimes; 
alphal=K (State, 1 ) ; 
alpha2=K (State,  2 )  ; 

alpha=0;  %  parameter  not  needed;  dummy  value  only 
beta=0;  %  parameter  not  needed;  dummy  value  only 
Aparam=0;  %  parameter  not  needed;  dummy  value  only 
Bparam=0;  %  parameter  not  needed;  dummy  value  only 
[k,T,To]  = 

phasetypeapproxmixed (A, alpha, beta, alphal , alpha2 , Aparam, Bparam, distribution) ; 

%  creates  phase-type  approximations  to  each  of  the  holding  time  distributions 


kl=k; 

T1=T; 

Tol=To; 

ze rovect or 1  =  zeros ( 1 ,  kl )  ; 

Psimatrixll  =  [  T1  Tol;  %  subgenerator  matrix  for  state  1 

zerovectorl  Qhat ( 1 , 1 ) *M] ; 


elseif  State  ==  2 

distribution  =  2; 

A  =  State2times; 
alphal=0;  %  parameter 
alpha2=0;  %  parameter 
alpha=K (State, 2) ; 
beta=K (State , 1 ) ; 

Aparam=0;  %  parameter 
Bparam=0;  %  parameter 
[k,T,To]  = 

phasetypeapproxmixed (A, alpha, beta, alphal , alpha2 , Aparam, Bparam, distribution) ; 

%  creates  phase-type  approximations  to  each  of  the  holding  time  distributions 


not 

not 


not 

not 


needed; 

needed; 


needed; 

needed; 


dummy 

dummy 


dummy 

dummy 


value 

value 


value 

value 


only 

only 


only 

only 


k2=k; 

T2=T; 

To2=To; 

ze rovect or 2  =  zeros ( 1 ,  k2 )  ; 

Psimatrix22  =  [  T2  To2;  %  subgenerator  matrix  for  state  2 

zerovector2  Qhat ( 2 , 2 ) *M] ; 


elseif  State  ==  3 

distribution  =  3; 

A  =  State3times; 

alpha=0;  %  parameter  not  needed;  dummy  value  only 
beta=0;  %  parameter  not  needed;  dummy  value  only 
alphal=K (State, 1 ) ; 
alpha2=K (State , 2 ) ; 

Aparam=0;  %  parameter  not  needed;  dummy  value  only 
Bparam=0;  %  parameter  not  needed;  dummy  value  only 
[k,T,To]  = 

phasetypeapproxmixed (A, alpha, beta, alphal , alpha2 , Aparam,  Bparam,  distribution)  ; 
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%  creates  phase-type  approximations  to  each  of  the  holding  time  distributions 
k3=k; 

T  3=T  ; 

To3=To; 

zerovector3=zeros ( 1 , k3 ) ; 

Psimatrix33  =  [  T3  To3;  %  subgenerator  matrix  for  state  3 

zerovector3  Qhat ( 3 , 3 ) *M] ; 


elseif  State  ==  4 

distribution  =  2; 

A  =  State4times; 
alpha=K (State, 2) ; 
beta=K (State ,  1 )  ; 

alphal=0;  %  parameter  not  needed;  dummy  value  only 
alpha2=0;  %  parameter  not  needed;  dummy  value  only 
Aparam=0;  %  parameter  not  needed;  dummy  value  only 
Bparam=0;  %  parameter  not  needed;  dummy  value  only 
[k,T,To]  = 

phasetypeapproxmixed (A, alpha, beta, alpha 1 , alpha2 , Aparam, Bparam, distribution) ; 

%  creates  phase-type  approximations  to  each  of  the  holding  time  distributions 


k4=k; 

T4=T  ; 

To4=To; 

zerovector4  =  zeros ( 1 ,  k4  )  ; 

Psimatrix44  =  [  T4  To4;  %  subgenerator  matrix  for  state  4 

zerovector4  Qhat ( 4 , 4 ) *M] ; 


elseif  State  ==  5 

distribution  =  5; 

A  =  State5times; 

Aparam=K (State, 1 ) ; 

Bparam=K (State, 2 ) ; 

alpha=0;  %  parameter  not  needed;  dummy  value  only 

beta=0;  %  parameter  not  needed;  dummy  value  only 

alphal=0;  %  parameter  not  needed;  dummy  value  only 

alpha2=0;  %  parameter  not  needed;  dummy  value  only 

[k,T,To]  = 

phasetypeapproxmixed (A, alpha, beta, alpha 1 , alpha2 , Aparam, Bparam, distribution) ; 

%  creates  phase-type  approximations  to  each  of  the  holding  time  distributions 


k5=k; 

T5=T; 

To5=To; 

zerovector5  =  zeros ( 1 ,  k 5 )  ; 

Psimatrix55  =  [  T5  To5;  %  subgenerator  matrix  for  state  5 

zerovector5  Qhat (5,  5) *M]  ; 


end 


end 


%  OFF-DIAGONAL  subgenerator  matrices  (see  p.3-15  of  thesis) 

%  ****IF  ADDITIONAL  STATES  ARE  ADDED,  THEY  MUST  BE  HARDCODED  HERE**** 

Psimatrixl2  =  zeros (kl+1 , k2+l) ; 

Psimat rix!2 (kl+1 , 1 ) =Qhat ( 1 , 2 ) *M; 


Psimatrixl3  =  zeros (kl+1 , k3+l) ; 
Psimat rix 13 (kl+1 , 1 ) =Qhat ( 1 , 3 ) *M; 

Psimatrixl4  =  zeros (kl+1 ,k4+l); 
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Psimatrixl4 (kl+1, 1) =Qhat (1, 4) *M; 

o. 

Psimatrixl5  =  zeros (kl+1 , k5+l ) ; 
Psimatrixl5 (kl+1, 1) =Qhat (1, 5) *M; 

o. 

o 

Psimatrix21  =  zeros (k2+l , kl+1 ) ; 
Psimatrix21 (k2+l , 1 ) =Qhat (2, 1 ) *M; 

o. 

o 

Psimatrix23  =  zeros (k2+l , k3+l ) ; 
Psimat r ix2  3 (k2  +  l , 1 ) =Qhat ( 2 , 3 ) *M; 

o. 

o 

Psimatrix24  =  zeros (k2+l , k4+l ) ; 
Psimatrix24 (k2+l , 1 ) =Qhat (2 , 4 ) *M; 

o. 

Psimatrix25  =  zeros (k2+l , k5+l ) ; 
Psimat rix2  5 (k2  +  l , 1 ) =Qhat ( 2 , 5 ) *M; 

g. 

o 

Psimatrix31  =  zeros (k3+l , kl+1 ) ; 
Psimatrix31 (k3+l , 1 ) =Qhat ( 3 , 1 ) *M; 

g. 

o 

Psimatrix32  =  zeros (k3+l , k2  +  l)  ; 
Psimatrix32 (k3+l , 1 ) =Qhat ( 3 , 2 ) *M; 

g. 

Psimatrix34  =  zeros (k3+l , k4+l ) ; 
Psimatrix34 (k3+l , 1 ) =Qhat ( 3 , 4 ) *M; 

g. 

o 

Psimatrix35  =  zeros (k3+l , k5+l)  ; 
Psimatrix35 (k3+l , 1 ) =Qhat ( 3 , 5 ) *M; 

g. 

o 

Psimatrix41  =  zeros (k4+l , kl+1 ) ; 
Psimatrix41 (k4+l, 1) =Qhat (4, 1) *M; 

g. 

o 

Psimatrix42  =  zeros (k4+l , k2+l ) ; 
Psimatrix42 (k4+l , 1 ) =Qhat (4,2) *M; 


Psimatrix43  =  zeros(k4+l,k3+l); 
Psimatrix43 (k4+l, 1) =Qhat (4, 3) *M; 

g. 

o 

Psimatrix45  =  zeros(k4+l,k5+l); 
Psimatrix45 (k4+l, 1) =Qhat (4, 5) *M; 

g. 

o 

Psimatrix51  =  zeros(k5+l,kl+l); 
Psimatrix51 (k5+l , 1 ) =Qhat ( 5 , 1 ) *M; 

g. 

o 

Psimatrix52  =  zeros(k5+l,k2+l); 
Psimatrix52 (k5+l , 1 ) =Qhat ( 5 , 2 ) *M; 

g. 

o 

Psimatrix53  =  zeros (k5+l , k3+l)  ; 
Psimatrix53 (k5+l , 1 ) =Qhat ( 5 , 3 ) *M; 

g. 

o 

Psimatrix54  =  zeros(k5+l,k4+l); 
Psimatrix54 (k5+l , 1 ) =Qhat ( 5 , 4 ) *M; 

g. 

o 

Psimatrix  =  [Psimatrixll  Psimatrixl2 
Psimatrix21  Psimatrix22 
Psimatrix31  Psimatrix32 
Psimatrix41  Psimatrix42 


Psimatrixl3 

Psimatrix23 

Psimatrix33 

Psimatrix43 


Psimatrixl4 

Psimatrix24 

Psimatrix34 

Psimatrix44 


Psimatrixl5 

Psimatrix25 

Psimatrix35 

Psimatrix45 
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Psimatrix51  Psimatrix52  Psimatrix53  Psimatrix54  Psimatrix55 ] ; 

o, 

%  This  section  creates  the  expanded  degradation  rate  matrix  Lambda 

%  ****IF  ADDITIONAL  STATES  ARE  ADDED,  THEY  MUST  BE  HARDCODED  HERE**** 

O, 

V=diag (R) ; 

Lll  =  diag (diag (V ( 1 , 1 ) *ones (kl+1 ) ) ) ; 

L12  =  zeros (kl+1 , k2+l ) ; 

L13  =  zeros (kl+1 , k3+l ) ; 

L14  =  zeros (kl+1 , k4+l ) ; 

L15  =  zeros (kl+1 , k5+l ) ; 

L21  =  zeros (k2+l , kl+1 ) ; 

L22  =  diag (diag (V (2, 2) *ones (k2  +  l) ) )  ; 

L23  =  zeros (k2+l , k3+l ) ; 

L24  =  zeros (k2+l , k4+l ) ; 

L25  =  zeros (k2+l , k5+l ) ; 

L31  =  zeros (k3+l , kl+1 ) ; 

L32  =  zeros (k3+l , k2+l ) ; 

L33  =  diag (diag (V (3, 3) *ones (k3+l) ) ) ; 

L34  =  zeros (k3+l , k4+l ) ; 

L35  =  zeros (k3+l , k5+l ) ; 

L41  =  zeros (k4+l , kl+1 ) ; 

L42  =  zeros (k4+l , k2+l ) ; 

L43  =  zeros (k4+l , k3+l ) ; 

L44  =  diag (diag (V ( 4 , 4 ) *ones (k4+l ) ) ) ; 

L45  =  zeros (k4+l , k5+l ) ; 

L51  =  zeros (k5+l , kl+1 ) ; 

L52  =  zeros (k5+l , k2+l ) ; 

L53  =  zeros (k5+l , k3+l ) ; 

L54  =  zeros (k5+l , k4+l ) ; 

L55  =  diag (diag (V (5, 5) *ones  (k5  +  l) ) ) ; 

o, 

L  =  [Lll  L12  L13  L14  L15;  %  expanded  degradation  rate  matrix 

L21  L22  L23  L24  L25; 

L31  L32  L33  L34  L35; 

L41  L42  L43  L44  L45; 

L51  L52  L53  L54  L 5 5 ] ; 

o, 

%  Display  the  results  of  the  phase-type  approximations 

o, 

fprintf  (  '  Size  of  generator  matrix  (Psimatrix)  :  %g  x  %g  \n',  length (Psimatrix) , 
length (Psimatrix) ) 
fprintf ( ' \n ' ) 

fprintf  (  '  Size  of  expanded  degradation  rate  matrix  (Lambda)  :  %g  x  %g  \n', 
length  (L),  length  (L) ) 
fprintf ( ' \n ' ) 

fprintf ( '  ******************************  \n  '  ) 
fprintf ( ' \n ' ) 

o, 

%  This  section  computes  and  plots  the  failure  time  distribution  of  the 
%  system  exposed  to  the  SMP  environment 

o, 

for  w  =  l:length(t) 

FailureTimeProb (w)  =  invt_lap_Solo (t (w) , Psimatrix,  L) ; 
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if  FailureTimeProb (w)  >  1 

FailureTimeProb (w)  =1;  %  eliminates  CDF  values  >  1 

elseif  FailureTimeProb (w)  <  0 

FailureTimeProb (w)  =0;  %  eliminates  CDF  values  <  0 

end 

end 

%  Compute  first  and  second  moment  of  failure  time  distribution 

o, 

First_Moment  =  invt_lap_f irst_moment (WearThreshold, Psimatrix, L) ; 

o, 

Second_Moment  =  invt_lap_second_moment (WearThreshold, Psimatrix, L) ; 

o, 

o, 

%  WRITE  OUTPUTS  TO  A  TAB  DELIMITED  TEXT  FILE 

O, 

fid  =  f open ( ' f : \AFIT\Thesis\Solo  Thesis\MATLAB  filesXSMP 
simulationXresult s  Xchemanalytical . out '  ,  ' w '  )  ; 
for  k=l : length (t) 

fprintf (fid,  '%4.6f\t%4.6f\n',t (k) , FailureTimeProb (k) ) ; 

end 

f close (fid) ; 
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Appendix  D.  Shared  Code  for  Examples 


%  Plot_SMP_f ive .m 

o, 

%  This  MATLAB  program  runs  the  semi-Markov  process  (SMP)  simulation 
%  ( semi_Markov_f ive . m)  and  then  the  'simulation  plus  phase-type  approximations 
%  plus  analytical  solution'  program  SMPfivestates .m.  The  CDF  values  are  then 
%  compared  in  a  plot  (Figure  2) . 

o, 

O, 

%  Author:  Captain  Christopher  J.  Solo,  M.S.  candidate,  Operations  Research 
%  Air  Force  Institute  of  Technology 

%  Date:  January  29,  2004 

o, 

o, 

clear  all; 


t=0 .05:0 .005 : 10; 

o, 

semi_Markov_f ive; 
obtain 

o, 

o, 

SMPfivestates ; 


%  vector  of  time  points  at  which  to  evaluate  both  CDFs 

%  calls  program  to  simulate  semi-Markov  process  (SMP)  and 
failure  times 

%  calls  program  to  simulate  SMP  (long  term),  compute  phase-type 
%  approximations,  &  produce  analytical  results  of  failure  time 
%  distribution 


%  F  =  cdf  values  from  simulation 

%  FailureTimeProb  =  cdf  values  from  phase-type  analytical  solution 

o, 

%  Produces  plot  of  simulated  failure-time  distribution  vs.  analytical 
%  CDF  derived  via  conversion  of  environment  process  from  SMP  to  CTMC 


figure  (2) ; 
plot (t , F, ' r ' ) ; 
hold  on; 

plot (t ,  FailureTimeProb,  ' b ' ) ; 
grid  off; 

title (' Failure-time  CDF  values:  Simulated  vs.  Analytical'); 
xlabel ( ' t '  )  ; 
ylabel ( 'P (T  <  t)  '  )  ; 

legend (' Simulated  CDF ',' Analytical  CDF',0); 
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PROGRAM  semi_Markov_f ive . m 

The  purpose  of  this  MATLAB  program  is  to  simulate  a  finite-state  semi-Markov 
process.  The  process  is  simulated  in  order  to  validate  analytical  solutions 
for  the  distribution  and  moments  of  the  failure  time  of  systems  with 
degradation  threshold  L  subject  to  a  SMP  environment.  The  program  uses 
function  "rando"  in  order  to  select  the  next  state  after  a  state  transition. 


Orig  Author : 

Date : 
Revised  by: 

Last  Revised: 

■k'k'k'k-k'k'k-k'k'k'k-k'k-k 


Jeffrey  P.  Kharoufeh,  Ph.D.  candidate,  IE  &  OR, 
Penn  State  University 
January  15,  2001 

Captain  Christopher  Solo,  M.S.,  OR, 

Air  Force  Institute  of  Technology 
January  18,  2004 


VARIABLE  DEFINTIONS 


WearThreshold  =  degradation  threshold  (system  fails  after  accumulating  this 
amount  of  damage) 
k  =  an  index  variable 

★★★★★it************************************************************************* 


%  VECTOR/FUNCTION  DEFINITIONS 

O, 

%  T  =  Vector  of  failure  time  observations 

%  B  =  Vector  for  the  initial  probability  distribution  of  the  semi-Markov  process 
%  P  =  Probability  transition  matrix 

%  V  =  Vector  of  degradation  rates  (i.e.,  V(i)=  degrad,  rate  when  environment  is 
%  in  state  i) . 

%  Z  =  Vector  of  environment  states 

o, 

%  The  program  assumes  a  state  space  of  the  form  S={1,2, . . .K} 

o, 

%  Simulation  of  semi-Markov  process  (SMP) 


semi_Markov_input_f ive;  %  Obtain  initialization  parameters 

Z  =  [];  %  Initialize  Z. 


for  k  =  1 : length (Lifetime) 

Z  =  El; 

Z ( 1 ) =rando (B) ; 


if  Z  ( 1 )  ==  S  (1 )  I  Z  ( 1 )  ==  S  (3) 

Lifetime(k)  =  betarnd (K ( Z ( 1 ) , 1 ) , K ( Z ( 1 ) , 2 ) ) ;  %  Time  spent  in  initial  state 
elseif  Z ( 1 )  ==  S (2 )  I  Z  (1 )  ==  S  (4 ) 

Lifetime(k)  =  weibrnd (K ( Z ( 1 ) , 1 ) , K ( Z ( 1 ) , 2 ) )  ; 
elseif  Z (1)  ==  S (5) 

Lifetime(k)  =  gamrnd (K ( Z ( 1 ) , 1 ) , K ( Z ( 1 ) , 2 ) ) ; 

end 

D  =  0.0;  %  At  time  0,  distance  traveled  is  0 

D  =  D  +  Lifetime (k) *V(Z (1) ) ;  %  Add  traveled  distance  to  cumulative  distance 

i  =  l; 

while  D  <  WearThreshold  %  Do  while  distance  traveled  is  less  than  WearThr 
Z(i+1)  =  rando (P (Z ( i ),:)) ;  %  Use  the  matrix  P  to  determine  next  state 

if  Z  (i+1)  ==  S ( 1 )  I  Z  (i  +  1)  ==  S  (3) 
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new_time  =  betarnd (K (Z (i  +  1) , 1) , K (Z (i  +  1 ) , 2 ) )  ; 

%  Generate  beta  rv  based  on  current  state 

elseif  Z  (i+1)  ==  S(2)  |  Z (i  +  1)  ==  S(4) 

new_time  =  weibrnd (K (Z (i  +  1) , 1) , K  (Z (i  +  1 ) , 2 ) ) ; 

%  Generate  beta  rv  based  on  current  state 
elseif  Z (i+1)  ==  S  (5 ) 

new_time  =  gamrnd (K ( Z ( i+1 ) , 1 ) , K ( Z (i+1 ) , 2 ) ) ; 

%  Generate  gamma  rv  based  on  current  state 
end 

D  =  D  +  new_time*V (Z (i+1) ) ;  %  Update  cumulative  damage 

Lifetime (k)  =  Lifetime (k)  +  new_time;  %  Update  total  lifetime 
i=i  +  l 

end 

Lif et ime (k) =Lifetime (k) - ( 1/V ( Z ( i ) ) ) * (D-WearThreshold) ;  %  Subtract  extra  time 

%  after  WearThr  is 
%  reached 

end 

get_cdf_f ive;  %  computes  the  empirical  distribution  of  observed  failure  times 
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phasetypeapproxmixed . m 

The  purpose  of  this  MATLAB  program  is  to  compare  the  values  of  an  arbitrary 
probability  distribution  with  those  of  its  2-phase  Coxian  or  k-phase 
generalized  Erlang  approximation. 


Author : 

Date : 
Last  Revised: 
References : 


Captain  Christopher  J.  Solo,  M.S.  candidate, 

Operations  Research,  Air  Force  Institute  of  Technology 
September  13,  2003 
January  18,  2004 

Perros,  H.G.  (1994).  Queueing  Networks  with  Blocking,  18-21. 
Altiok,  T.  (1985) .  On  the  phase-type  approximations  of  general 
distributions,  HE  Transactions,  Vol.  17,  No.  2,  110-116 
Neuts,  M.F.  (1981)  .  Matrix-Geometric  Solutions  in  Stochastic 
Models,  41-63. 

Sauer,  C.H.  and  Chandy,  K.M.  (1975)  .  Approximate  analysis  of 
central  server  models,  IBM  Journal  of  Research  and  Development, 
Vol.  19,  No.  3,  301-313. 

Cox,  S.R.  (2003) .  Draft  MATLAB  code,  Air  Force  Institute  of 
Technology . 

Law,  A.M.  and  Kelton,  W.D.  (2000)  .  Simulation  Modeling  and 
Analysis . 

Wackerly,  D.,  Mendenhall,  W.,  and  Scheaffer,  R.  (2002) . 
Mathematical  Statistics  with  Applications,  444. 


function  [k,T,To]  = 

phasetypeapproxmixed (A, alpha, beta, alphal , alpha2 , Aparam, Bparam, distribution) 


CDFdistr=[] ; 
pdfdistr= [ ] ; 
CDFapprox= [ ] ; 
pdf approx= [ ] ; 


CDF  values  of  probability  distribution 
density  values  of  probability  distribution 
CDF  values  of  2-phase  Coxian  approximation 
density  values  of  2-phase  Coxian  approximation 


INDICATE  TYPE  OF  INPUT  DISTRIBTUION  (if  known) 


1  =  exponential 

2  =  Weibull 

3  =  beta 

4  =  uniform 

5  =  gamma 


syms  x;  %  symbolic  x  used  in  calculation  of  first  three  moments 
%  INDICATE  INPUT  DISTRIBUTION 

%  EXPONENTIAL  DISTRIBUTION 


if  distribution  ==  1 

lambda  =  27.74563544200000;  %  parameter  of  exponential  distribution 

f  =  lambda*exp (-lambda*x) ;  %  exponential  pdf 

obs  =  2000;  %  number  of  observations  from  distribution 

t  =  linspace ( 0 , . 5 , 5*obs ) ; 

%  time  points — 'obs'  points  between  0  and  2,  including  0 


D-4 


%A  =  RANDOM (' exp 1/lambda,  obs,  1 ) 
observations 


generates  a  'obs'  x  1  vector  A  of 


%  WEIBULL  DISTRIBUTION 

O, 

elseif  distribution  ==  2 

f  =  alpha* (betaA (-alpha) )* (xA (alpha  -  1) ) *exp (- (x/beta) Aalpha) ;  %  Weibull  pdf 
obs  =  200; 

t  =  linspace  ( 0 , 1 . 4 , 5*obs ) ;  %  time  points-- ' obs '  points  between  0  and  2 
%A  =  RANDOM (' weib ', beta, alpha, obs , 1 ) ;  %  generates  'obs'  x  1  vector  of  obs 

o, 

9^'k-k-k-k'k-k'k'k-k'k-k-k'k-k-k-k-k-k'k-k-k-k'k-k'k-k-k-k-k-k'k-k-k-k-k-k'k-k-k-k-k-k'k 

%  BETA  DISTRIBUTION 

O, 

elseif  distribution  ==  3 

BetaParam  =  exp (gamma In (alphal ) +gammaln ( alpha 2 ) -gammaln (alphal+alpha2 ) )  ; 
f  =  ( (xA (alphal  -  1))*(1  -  x) A (alpha2  -  1) ) /BetaParam; 

%  beta  pdf  from  Law  and  Kelton 
obs  =  200; 

t  =  linspace  (0, 1, 5*obs) ; 

%  time  points — 'obs'  points  between  0  and  2,  including  0 
%A  =  RANDOM ( 'beta ', alphal, alpha2, 5 *obs, 1) ; 

%  generates  a  'obs'  x  1  vector  A  of  observations 

o, 

9^'k-k-k-k'k-k'k'k-k'k-k-k-k-k-k-k-k-k'k-k-k-k'k-k'k-k-k-k-k-k'k-k-k-k-k'k'k-k-k-k-k-k'k 

%  UNIFORM  DISTRBUTION 

O, 

elseif  distribution  ==  4 

f  =  1/  (theta2  -  thetal) ; 
obs  =  2000; 

%A  =  RANDOM (' unif ', thetal , theta2 , 5*obs, 1 ) ;  %  generates  a  'obs'  x  1  vector  A 

%  of  observations 

O, 

9^-k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  ~k  -k  -k  -k  -k  -k  -k  -k  -k 

%  GAMMA  DISTRIBUTION 

O, 

elseif  distribution  ==  5 
obs  =  2000; 

t  =  linspace ( 0 , 1 . 4 , 5*obs )  ; 

end 

%  %  Moments  of  given  distribution 

%  (use  this  when  input  distribution  is  known) 


o  if  (distribution  ==*  1)  |  (distribution 

’s  lowerlimit  =  0; 

o  upperlimit  =  inf; 

b  t  =  linspace ( 0 , 2 , 500 )  ; 

Including  0 

o  elseif  distribution  ==  3 
b  lowerlimit  =  0; 

b  upperlimit  =  1; 

b  t  =  linspace ( 0 , 2 , 500 )  ; 

b  elseif  distribution  ==  4 
b  lowerlimit  =  thetal; 

b  upperlimit  =  theta2; 


2) 


exponential  or  Weibull 


%  upper  limit  of  integration 
time  points--500  points  between  0  and  2, 

%  beta  distribution 
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%  t  =  linspace ( 0 , theta2  +  4,500); 

%  end 
2-  2- 

%  ml  =  eval ( int ( (x* f) , lowerlimit , upper limit )) ;  %  1st  moment  of  given 

%  distribution 

%  m2  =  eval ( int (( (xA2 )* f) , lowerlimit , upperlimit )) ;  %  2nd  moment  of  given 

%  distribution 

%  m3  =  eval ( int (( (xA3 )* f) , lowerlimit , upperlimit )) ;  %  3rd  moment  of  given 

%  distribution 


%  Moments  derived  from  input  data 

o, 

%if  A  is  the  vector  of  input  values  (holding  time  observations) 
for  m  =  1 : length (A) 

cubed(m,l)  =  (A(m,l))A3; 

end 


ml  =  mean (A) ; 

m2  =  var (A)  +  mlA2; 

m3  =  ( 1/length (A) )* (sum (cubed) ) ;  %  from  Wackerly 


%  Moments  derived  via  successive  derivatives  of  Laplace  transform 
%  THIS  APPLIES  ONLY  TO  THE  EXPONENTIAL  DISTRIBUTION 

O, 

%  ml=l/lambda; 

%  m2=2 /lambda A2 ; 

%  m3=6/lambdaA3 ; 

o, 

%  (Squared)  Coefficient  of  Variation 

o, 

c  =  (m2-ml A2 ) /ml A2 ;  %  squared  coefficient  of  variation  =  var/meanA2 
z  =  3* (m2  A2 ) ; 
v  =  2*ml*m3; 

o, 

S-icici^ic^ificifif^'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

%  FOR  c  >  1 

%  approximation  is  made  using  the  2-phase  Coxian  distribution 

o, 

if  c  >  1 

k=2;  %  default  number  of  phases 
phase  =  '3-moment,  2-phase  Coxian'; 

%  Moments  of  2-phase  Coxian  approximation 

o, 

Y  =  (6*ml  -  (3*m2/ml) ) / ( (6*m2A2/4*ml)  -  m3);  %  Altiok  (1985),  p.112 
X  =  (1/ml)  +  ( (m2*Y) / (2*ml) ) ;  %  Perros  (1994),  p.27 

mul  =  (X  +  sqrt (XA2  -  4*Y))/2;  %  first  moment  of  2-phase  approximation 

mu2  =  X  -  mul;  %  second  moment  of  2-phase  approximation 
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Sr-k-k-k-k-k-k-k-kk'k'k'k-k'k'k'k-k-k'k'k'k-k'k'k'k'k-k'k'k'k-k-k'k'k'k-k'k'k'k'k-k'k'k'k-k-k'k'k'k-k'k'k'k'k-k'k'k'k-k-k'k'k'k-k'k'k'k'k-k'k'k'k-k-k 


%  Underlying  CTMC 

o, 

a  =  (mu2/mul ) * ( (ml *mul )  -  1);  %  prob  of  transition  from  phase  1  to  phase  2 
T  =  [-mul  a*mul;  %  matrix  of  rate  transitions  among  2  phases 

0  -mu2 ] ; 

To  =  [(l-a)*mul;  mu2] ;  %  2x1  vector  of  rates  of  transition  out  of  both  states 

%  into  absorbing  state 

initialprob  =[10];  %  initial  probability  vector 

ones=[l;l];  %  vector  of  ones 


%  Original  distribution  and  approximation  values  and  plots 


PHvalues (c, k, T, To, phase, t, initialprob, ones, distribution, alpha, beta, alphal, alpha2 
, Aparam, Bparam) ;  %  calls  PHvalues. m  file 

o, 

9^-k-k'k'k-k-k-k-k-k-k-k-k-k'k-k-k'k'k-k-k-k'k'k-k-k-k'k-k-k-k-k'k'k-k-k-k-k-k-k-k-k'k-k-k'k'k-k-k-k-k'k-k-k-k-k-k-k-k'k'k'k-k-k-k-k-k-k-k'k'k-k-k'k'k 

%  FOR  0.5  <=  c  <=  1 

%  use  the  two-moment  approximation  from  above  (see  Perros  (1994),  p.30) 


elseif  (c  >=  0.5)  &  (c  <=  1) 
k=2;  %  default  number  of  phases 
mul  =  2/ml; 
mu2  =  1/ (ml*c) ; 

phase  =  '2-moment  2-phase  Coxian'; 

o, 

%  Underlying  CTMC 


a  =  1/  (2*c) ; 

T  =  [-mul  a*mul; 

0  -mu2 ] ; 

To  =  [ (1-a) *mul;  mu2]; 
into  absorbing  state 
initialprob  =  [10]; 
ones  =  [ 1 ;  1 ] ; 


probability  of  transition  from  phase  1  to  phase  2 
matrix  of  rate  transitions  among  2  phases 

i  2x1  vector  of  rates  of  transition  out  of  both  states 

%  initial  probability  vector 
i  vector  of  ones 


Sfe  "k  "k  -k  -k  k  ★  "k  -k  k  ★  ★  k  *  k  ★  "k  -k  -k  k  "k  k  -k  k  "k  ~k  -k  -k  k  "k  -k  -k  -k  k  "k  it  -k  k  ★  ★  "k  -k  k  ★  ★  -k  -k  k  "k  "k  -k  k  "k  ★  ★  -k  k  "k  "k  ★  -k  k  ★  ★  -k  k  ★ 

%  Original  distribution  and  approximation  values  and  plots 


PHvalues (c,  k,  T, To,  phase, t , initialprob, ones, distribution, alpha, beta, alphal , alpha2 
, Aparam, Bparam) ;  %  calls  PHvalues. m  file 

o, 

^■k-k-k-k-k'k-k-k-k-k'k-k'k-k'k-k-k-k-k'k-k'k-k'k'k-k-k-k'k-k-k-k-k'k-k-k-k'k'k-k-k-k'k-k-k-k-k'k-k-k-k'k'k-k-k-k'k-k-k'k-k'k-k-k-k'k 
^  ■ k'k-k-k'k'k'k-k'k'k'k'k'k'k'k'k-k-k'k'k'k-k'k'k'k'k-k'k'k'k-k-k'k'k'k-k'k'k'k'k-k'k'k'k-k-k'k'k'k-k'k'k'k'k'k'k'k'k-k-k'k'k'k-k'k'k 

%  FOR  c  <  0 . 5 

%  approximation  is  made  using  a  generalized  Erlang  distribution  with  k  phases 

%  (see  Perros  (1994),  p.29-30) 

else 

o,  o, 

%  Linear  program  to  determine  minimum  number  of  phases  in  generalized  Erlang 
%  distribution 

o, 

Aeq= [ ] ; 
beq= [ ] ; 
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ub=  [  ] ; 

d  =  [1];  %  choose  k  such  that  1/k  <=  c  AND  l/(k-l)  >=  c 

Y  =  [-  :  ; 

l]  ; 

b  =  [-1/c  (1/c) +1] ; 
lb  =  0; 

k  =  ceil ( linprog (d, Y, b, Aeq, beq, lb, ub) ) ;  %  chooses  first  integer  greater  than 
solution 

phase  =  'k-phase  generalized  Erlang'; 

o, 

%  Underlying  CTMC 

o, 

%  prob  of  transition  from  phase  i  to  phase  i+1  (see  Ferros  (1994),  p.29-30) 
a  =  1  -  ( ( (2*k*c) +k-2- (kA2+4- (4*k*c) ) A (1/2) ) / (2* (c+1) * (k-1) ) ) ; 

o, 

mu  =  (1+ ( (k-1) *a) ) /ml;  %  rate  of  transition 

o, 

3^~k  ic  ~k  ~k  i<  ~k  ~k  ~k  i<  ~k  •)(  ~k  ~k  i<  ~k  ~k  ic  ic  i<  ~k  is  ~k  i<  ~k  ~k  ~k  ic  i<  ~k  ~k  ~k  -k  i<  ~k  ~k  ~k  i<  •)<  ~k  ~k  ~k  ic  ~k  ~k  ic  ic  i<  ~k  ic  ~k  i(  ~k  ~k  ~k  ~k  i(  •)(  ~k  ic  ~k  i<  ~k  ~k  ~k  i<  ~k  ~k 
o, 

%  This  "for"  loop  creates  the  kxk  matrix  T  of  transition  rates 

o, 

T  =  zeros (k,k);  %  preallocates  zeroes  to  the  matrix 
T  ( 1 , 2 ) =a*mu; 
for  ( f =1 : k) 

for  (g=l:k) 


if  f  ==  g 

T  ( f , g )  =  -mu ; 
elseif  f  >  g 

T  (f,  g)  =  0; 
elseif  g  ==  (f  +  1) 
T  ( f ,  g)  =  mu;  % 
elseif  g  >  (f  +  1) 

T  (f,  g)  =  0; 

end 


%  diagonal  elements 
%  below  diagonal  elements 
just  above  diagonal  elements 
%  above  diagonal  elements 


end 

end 

T ( 1 , 2 ) =a*mu; 

To  =  zeros  (k,l);  %  pre-allocates  zeroes  to  the  vector 
ones  =  zeros (k,l); 
initialprob  =  zeros (l,k); 

To  (1,1)  =  (l-a)*mu;  %  can  enter  absorbing  state  from  first  phase  with  prob  (1-a) 
To(k,l)  =  mu;  %  rate  of  transition  out  of  last  phase  into  absorbing  state 
ones (1,1) =1 ; 
ones (k, 1 ) =1 ; 
for  r=2 : k-1 

To(r,l)  =  0;  %  kxl  vector  of  rates  of  transition  out  of  phases  2  through  k- 

%  1  into  absorbing  state 

ones  (r, 1)  =1;  %  kxl  vector  of  ones  (1,1,1,  .  .  .1) 

end 

initialprob ( 1 , 1 )  =  1;  %  initial  probability  vector  (1,0,0,  .  .  .0) 

for  s=2:k 

initialprob (1, s)  =  0; 

end 
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o\°  o\°  >  hd  o\°  o\°  o\° 


Original  distribution  and  approximation  values  and  plots 


Hvalues (c, k, T, To, phase, t 
Aparam, Bparam) ;  %  calls 


,  initialprob, ones, distribut 
PHvalues.m  file 


ion, alpha, beta, alphal 


alpha2 


★  ★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★■St* 

end 
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PHvalues  .m 


The  purpose  of  this  MATLAB  program  is  to  compute  the  approximated 
CDF  values,  the  exact  CDF  values,  and  their  differences,  as  derived 
from  the  phasetypeapproxmixed . m  program.  This  program  then  plots  the 
approximated  CDF  vs .  exact  CDF  values . 


Author : 

Date  : 
Last  Revised: 
References : 


Captain  Christopher  J.  Solo,  M.S.  candidate. 

Operations  Research,  Air  Force  Institute  of  Technology 
September  13,  2003 
January  18,  2004 

Perros,  H.G.  (1994).  Queueing  Networks  with  Blocking,  18-21. 
Neuts,  M.F.  (1981)  .  Matrix-Geometric  Solutions  in  Stochastic 
Models,  41-63. 


function  []  = 

PHvalues (c, k, T, To, phase, t, initialprob, ones, distribution, alpha, beta, alphal, alpha2 
, Aparam, Bparam) 

%  Output  Display 

o, 

fprintf ( ' \n ' ) 

%fprintf  (  '  State:  %g  \n',  State) 
fprintf ( ' \n ' ) 

fprintf  (  '  Number  of  phases:  %g  \n',k) 
fprintf ( ' \n ' ) 

fprintf  (  '  Type  of  approximation:  %s  \n',  phase) 
fprintf ( ' \n 1 ) 

fprintf  (  '  Coefficient  of  variation  c  =  %f  \n',  c) 
fprintf ( ' \n ' ) 

fprintf  (  '  ******************************  \n ' ) 
fprintf ( ' \n ' ) 


%  fprintf ('  t  approx  CDF  exact  CDF  Difference  \n') 

o. 


for  u=l : length (t ) 

if  distribution  ==  1 

CDFdistr(u)  =  expcdf (t (u) , 1 /lambda) ; 
pdfdistr(u)  =  exppdf (t (u) , 1 /lambda) ; 
trueCDF(u,l)  =  t (u)  ; 
trueCDF(u,2)  =  CDFdistr (u) ; 
elseif  distribution  ==  2 

CDFdistr (u)  =  weibcdf (t (u) , beta, alpha) ; 

pdfdistr(u)  =  weibpdf (t (u)  ,  beta,  alpha)  ; 
trueCDF(u,l)  =  t (u)  ; 
trueCDF(u,2)  =  CDFdistr (u) ; 
elseif  distribution  ==  3 

CDFdistr (u)  =  betacdf (t (u) , alphal , alpha2 ) ; 

pdfdistr(u)  =  betapdf (t (u) , alphal , alpha2 ) ; 
trueCDF(u,l)  =  t (u)  ; 
trueCDF(u,2)  =  CDFdistr (u) ; 
elseif  distribution  ==  4 

CDFdistr (u)  =  unif cdf (t (u) , thetal, theta2 )  ; 

pdfdistr(u)  =  unifpdf (t (u)  ,  thetal,  theta2 )  ; 


%  CDF  of  the 
%  pdf  of  the 


%  CDF  of  the 
%  pdf  of  the 


%  CDF  of  the 
%  pdf  of  the 


%  CDF  of  the 
%  pdf  of  the 


distribution 

distribution 


distribution 

distribution 


distribution 

distribution 


distribution 

distribution 


D-10 


trueCDF(u,l)  =  t (u) ; 
trueCDF(u,2)  =  CDFdistr (u) ; 
elseif  distribution  ==  5 

CDFdistr (u)  =  gamcdf (t (u) , Aparam, Bparam) ;  %  CDF  of  the  distribution 

pdfdistr(u)  =  gampdf (t (u) , Aparam, Bparam) ;  %  pdf  of  the  distribution 

trueCDF(u,l)  =  t (u) ; 
trueCDF(u,2)  =  CDFdistr (u)  ; 

end 

CDFapprox(u)  =  1-initialprob* (expm (T*t (u) ) ) ‘ones; 

%  CDF  of  approximation  (see  Neuts,  1981) 

pdfapprox(u)  =  initialprob* ( (expm (T* (t (u) ) ) ) *To) ; 

%  pdf  of  approximation  (see  Perros,  1994) 

%  fprintf('  %f  %f  %f  %f\n',  t (u) ,  CDFdistr (u) , 

CDFapprox  (u) ,  abs (CDFdistr (u) -CDFapprox (u) ) ) 
end 

9-ieicieieieieieiciiciiciciiciciicici?iicicicicicicicicic*^**'k-k'k'k'k-k'k-k'k-k-k-k'k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k'k'k'k'k'k 

%  Kolmogorov-Smirnov  test 

%  (determines  if  approximated  values  and  exact  values  come  from  same 
%  distribution) 

H  =  kstest2 (CDFapprox, CDFdistr )  ; 
value  =  zeros (1, length (t) ) ; 
for  j  =  1: length (t) 

value  (j)  =  (CDFapprox ( j ) -CDFdistr  (  j )) A2 ; 

end 

%res  =  0 . 5*sum (value) 

%H  %  H  =  1  =>  reject  null  that  both  come  from  same  distribution 

%c 

%k 

%phase 

^ieicieieieieieiciiciiciciiciiciiciciiciicicicicicicicicic^**'k'k'k'k'k'k-k-k-k'k-k-k-k'k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k'k'k'k'k'k 

%  Plot  of  k-phase  PH  approx  CDF  vs.  exact  CDF 

figure ( 1 ) ; 

plot  (t, CDFapprox,  '  r  :  '  )  ; 
hold  on; 

plot  (t, CDFdistr,  ' b ' )  ; 
grid  off; 

title  ('CDF  values:  Phase-type  vs.  Input'); 

xlabel  (  ' t '  )  ; 

ylabel ( ' P (T  <  t)  '  )  ; 

legend ( ' Phase-type ' , ' Input ' , 0 ) ; 

2-  -k'k'k'k'k'k'k-k-k-k-k-k-k-k-k-k’k-k'k'k'k-k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k-k'k-k-k-k'k-k-k-k'k'k-k'k-k'k-k'k-k'k'k'k'k'k'k 
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PROGRAM  invt_lap_Solo . m 

The  purpose  of  this  MATLAB  program  is  to  approximate  the  inverse  transform  of 
a  one-dimensional  Laplace  transform  in  order  to  find  the  moments  of  the 
probability  distribution,  G(t) .  The  program  is  based  on  the  algorithm  of  Abate 
and  Whitt  (1995)  . 


Orig  Author: 

Date : 
Revised  by: 

Date : 
References : 


Jeffrey  P.  Kharoufeh,  Ph.D.  candidate,  IE  &  OR, 

Penn  State  University 
January  23,  2001 

Captain  Chris  Solo,  M.S.  candidate,  OR, 

Air  Force  Institute  of  Technology 
January  29,  2004 

Abate,  J.  and  W.  Whitt  (1995)  .  Numerical  Inversion  of  the 

Laplace  Transform  of  Probability  Distribution.  ORSA  Journal  on 
Computing,  7,  36-43. 


%Initialize  variables,  set  parameters 

function  fl  =  invt_lap_Solo (t, Psimatrix, L)  %  inputs  time  vector,  generator 

%  matrix,  and  degradation  matrix 

rho=0.8;  qx= [ 0 . 8 ] ;  tx=[0];  m=ll;  c= [ ] ;  ga=8;  A=ga*log (10) ;  mm=2Am; 

O, 

for  k=0:m 

d=nchoosek  (m,  k)  ; 
c= [c  d]  ; 

end 

for  t  =  t; 
tx  =  t; 
ntr=15 ; 
u=exp (A/2 ) / t ; 
x= A/ (2*t) ; 
h=pi/t; 

su=zeros (m+2 ) ; 

sm=f ailLST (x,  0 ,  Psimatrix, L) / 2 ; 
for  k=l:ntr 
y=k*h; 

sm=sm+ ( (-1) Ak) *failLST (x,y,  Psimatrix,  L)  ; 

end 

su ( 1 ) =sm; 
for  k=l : 12 
n=ntr+k; 
y=n*h; 

su(k+l)=su(k)+( (-1) An) *failLST(x,y, Psimatrix, L) ; 

end 

avl=0;  av2=0; 
for  k=l : 12 

avl=avl  +  c (k) *su  (k) ; 
av2=av2  +  c (k) *su  (k+1) ; 

end 

fl  =  u*avl/mm;  f2=u*av2/mm;  qx=[qx  f2]; 

end 
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%  Orig  Author:  Jeffrey  P.  Kharoufeh,  Ph.D.  candidate,  IE  &  OR, 

%  Penn  State  University 

%  Date:  June  2000 

^•k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

function  eq=failLST (x, y, Psimatrix, L) 

%  This  program  computes  the  LST  of  the  failure  time  according  to 
%  Kharoufeh' s  paper,  Sipe's  results 

s  =  x+y*i; 

I=eye (size (Psimatrix) ) ;  %  Create  identity  matrix 

Al=inv (L) *( (Psimatrix- (s* I )) *5 . 0 ) ;  %  The  last  argument  is  the  degradation 

%  threshold  level,  x. 

A2=expm (A1 ) ; 

J  =  zeros (1, (length (Psimatrix) -1) ) ;  %  total  number  of  phases  minus  one  =  number 

%  of  zeros  to  use 

alphavec  =  [1  J] ;  %  initial  probability  vector:  Ensure  the  #  columns  =  # 

%  states  (phases) . 
m=ones (size (Psimatrix, 1) , 1) ; 
z  =  ( 1/s ) * (alphavec*A2 *m) ;  %  get  the  cdf 

eq  =  real(z);  %  get  the  cdf 
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PROGRAM  invt_lap_first_moment .m 

The  purpose  of  this  MATLAB  program  is  to  approximate  the  inverse  transform  of 
a  one-dimensional  Laplace  transform  in  order  to  find  the  moments  of  the 
probability  distribution,  G(t) .  The  program  is  based  on  the  algorithm  of  Abate 
and  Whitt  (1995)  . 


Orig  Author: 

Date : 
Revised  by: 

Date : 
References : 


Jeffrey  P.  Kharoufeh,  Ph.D.  candidate,  IE  &  OR, 

Penn  State  University 
January  23,  2001 

Captain  Chris  Solo,  M.S.  candidate,  OR, 

Air  Force  Institute  of  Technology 
January  29,  2004 

Abate,  J.  and  W.  Whitt  (1995)  .  Numerical  Inversion  of  the 

Laplace  Transform  of  Probability  Distribution.  ORSA  Journal  on 
Computing,  7,  36-43. 


Initialize  variables,  set  parameters 


function  fl  =  invt_lap_first_moment (t, Psimatrix, L)  %  inputs  time  vector, 

%  generator  matrix, 

%  and  degradation  matrix 

rho=0.8;  qx= [ 0 . 8 ] ;  tx=[0];  m=ll;  c= [ ] ;  ga=8;  A=ga*log (10) ;  mm=2Am; 

O, 

for  k=0:m 

d=nchoosek  (m,  k)  ; 
c= [c  d]  ; 

end 

for  t  =  t; 
tx  =  t; 
ntr=15 ; 
u=exp (A/2 ) / t ; 
x=A/ (2*t) ; 
h=pi/t; 

su=zeros (m+2 ) ; 

sm=f irst_momentLST (x, 0 , Psimatrix, L) /2 ; 
for  k=l : ntr 
y=k*h; 

sm=sm+ ( (-1 ) Ak) *f irst_momentLST (x, y, Psimatrix, L) ; 

end 

su ( 1 ) =sm; 
for  k=l : 12 
n=ntr+k; 
y=n*h; 

su(k+l)=su(k)  +  (  (-1)  An)  *f  irst__momentLST  (x,  y,  Psimatrix,  L)  ; 

end 

avl=0;  av2=0; 
for  k=l : 12 

avl=avl  +  c (k) *su  (k) ; 
av2=av2+c (k) *su (k+1) ; 

end 

fl  =  u*avl/mm;  f2=u*av2/mm;  qx=[qx  f2]; 

end 
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PROGRAM  invt_lap_second_moment .m 

The  purpose  of  this  MATLAB  program  is  to  approximate  the  inverse  transform  of 
a  one-dimensional  Laplace  transform  in  order  to  find  the  moments  of  the 
probability  distribution,  G(t) .  The  program  is  based  on  the  algorithm  of  Abate 
and  Whitt  (1995) . 


Orig  Author: 

Date : 
Revised  by: 

Date : 
References : 


Jeffrey  P.  Kharoufeh,  Ph.D.  candidate,  IE  &  OR, 

Penn  State  University 
January  23,  2001 

Captain  Chris  Solo,  M.S.  candidate,  OR, 

Air  Force  Institute  of  Technology 
January  29,  2004 

Abate,  J.  and  W.  Whitt  (1995)  .  Numerical  Inversion  of  the 

Laplace  Transform  of  Probability  Distribution.  ORSA  Journal  on 
Computing,  7,  36-43. 


Initialize  variables,  set  parameters 


function  fl  =  invt_lap_second_moment (t , Psimatrix, L)  %  inputs  time  vector, 

%  generator  matrix, 

%  and  degradation  matrix 

rho=0.8;  qx= [ 0 . 8 ] ;  tx=[0];  m=ll;  c= [ ] ;  ga=8;  A=ga*log (10) ;  mm=2Am; 

O, 

for  k=0:m 

d=nchoosek  (m,  k)  ; 
c= [c  d]  ; 

end 

for  t  =  t; 
tx  =  t; 
ntr=15 ; 
u=exp (A/2 ) / t ; 
x=A/ (2*t) ; 
h=pi/t; 

su=zeros (m+2 ) ; 

sm=second_momentLST (x,  0 ,  Psimatrix,  L)  / 2 ; 
for  k=l : ntr 
y=k*h; 

sm=sm+ ( (-1 ) Ak) *second_momentLST (x, y, Psimatrix, L) ; 

end 

su ( 1 ) =sm; 
for  k=l : 12 
n=ntr+k; 
y=n*h; 

su(k+l)=su(k)+( (-1) An) *second_momentLST (x, y, Psimatrix, L) ; 

end 

avl=0;  av2=0; 
for  k=l : 12 

avl=avl  +  c (k) *su  (k) ; 
av2=av2+c (k) *su (k+1) ; 

end 

fl  =  u*avl/mm;  f2=u*av2/mm;  qx=[qx  f2]; 

end 
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%  PROGRAM  f irst_momentLST . m 

O, 

%  This  program  computes  the  LST  of  the  first  moment  according  to 
%  Kharoufeh  (2003)  and  Kharoufeh  and  Sipe  (2004) . 

o, 

%  Orig  Author:  Jeffrey  P.  Kharoufeh,  Ph.D.  candidate,  IE  &  OR, 

%  Penn  State  University 

%  Revised  by:  Captain  Christopher  Solo,  M.S.  candidate,  OR, 

%  Air  Force  Institute  of  Technology 

%  Last  Revised:  March  11,  2004 

O, 

function  eq=f irst_momentLST (x, y, Psimatrix,  L) 

clear  r ; 
n  =  1; 
r  =  x+y*i; 

Bl= (r*L) -Psimatrix; 

J  =  zeros  (1,  (length (Psimatrix) -1) ) ;  %  total  number  of  phases  minus  one  =  number 

%  of  zeros  to  use 

alphavec  =  [1  J] ;  %  initial  probability  vector:  Ensure  the  #  columns  =  # 

%  states  (phases) . 
m=ones (size (Psimatrix, 1) , 1) ; 

z  =  (1/r) * (factorial (n) ) *alphavec* ( (Bl) A (-n) ) *m;  %  get  the  first  moment 

eq  =  real(z);  %  get  the  cdf 
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%  PROGRAM  second_momentLST . m 

O, 

%  This  program  computes  the  LST  of  the  second  moment  according  to 
%  Kharoufeh  (2003)  and  Kharoufeh  and  Sipe  (2004) . 

o, 

%  Orig  Author:  Jeffrey  P.  Kharoufeh,  Ph.D.  candidate,  IE  &  OR, 

%  Penn  State  University 

%  Revised  by:  Captain  Christopher  Solo,  M.S.  candidate,  OR, 

%  Air  Force  Institute  of  Technology 

%  Last  Revised:  March  11,  2004 

O, 

function  eq=second_momentLST (x,  y,  Psimatrix,  L) 

%  This  program  computes  the  LST  of  the  failure  time  according  to 
%  Kharoufeh 's  paper,  Sipe's  results 
clear  r ; 
n  =  2 ; 
r  =  x+y*i; 

Bl=  (r*L) -Psimatrix; 

J  =  zeros  (1,  (length (Psimatrix) -1) ) ;  %  total  number  of  phases  minus  one  =  number 

%  of  zeros  to  use 

alphavec  =  [1  J] ;  %  initial  probability  vector:  Ensure  the  #  columns  =  # 

%  states  (phases) . 
m=ones (size (Psimatrix, 1) , 1) ; 

z  =  (1/r) * (factorial (n) ) *alphavec* ( (Bl) A (-n) ) *m;  %  get  the  first  moment 

eq  =  real(z);  %  get  the  cdf 
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%  PROGRAM  rando.m 

O, 

%  This  program  generates  a  random  variable  in  1,2, . . . ,n  given  a  distribution 
%  vector. 

o, 

%  Orig  Author:  Jeffrey  P.  Kharoufeh,  Ph.D.  candidate,  IE  &  OR, 

%  Penn  State  University 

%  Date:  June  2000 

o, 

function  [index]  =  rando (p)  %  Example  x(l)  =  rando(mu)  then  p  =  mu 

u  =  rand; 
i  =  1; 

s  =  p(l);  %  s  =  first  probability  in  p 

while  ( (u  >  s)  &  (i  <  length(p))),  %  Random  draw  compared  to 

probability 
i=i+l ; 
s=s+p ( i ) ; 
end 

index=i ; 


Must  be  a  stochastic  vector  (sum  =  1) 


Returns  state  to  transition  to  next 
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